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ABSTRACT 


The  goal  of  this  study  is  to  investigate  the  nonlinear  dynamics  of  two  ships  in 
close  proximity  towing.  The  sway  and  yaw  dynamics  of  both  the  leading  and  the  trailing 
ships  were  included  in  the  formulation.  Previous  studies  were  restricted  to  a  linear 
analysis,  which  can  accurately  predict  the  regions  of  stability  and  instability  for  the 
system.  The  mechanism  of  loss  of  stability  can  be  assessed  with  a  systematic  nonlinear 
analysis.  The  analysis  is  based  on  Taylor  series  expansions  of  the  equations  of  motion  up 
to  third  order  terms.  It  is  shown  that  the  primary  loss  of  stability  occurs  in  the  fonn  of 
Hopf  bifurcations  to  periodic  solutions.  A  nonlinear  stability  coefficient  was  calculated 
which  allows  characterization  of  the  stability  properties  of  the  resulting  limit  cycles.  The 
results  indicate  the  effects  of  ship  separation  and  towing  tension  on  motion  stability. 
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I.  INTRODUCTION 


One  of  the  main  concerns  in  both  naval  architecture  and  maneuvering  control  is 
the  stability  of  the  ship  towing.  In  particular,  the  close  proximity  ship-towing  problem 
has  many  unexplored  areas  of  investigation  ranging  from  building  a  mathematical  model 
of  the  two  towed  ships,  analyzing  the  directional  stability  of  the  combined  system,  local 
stability  of  the  maneuvering,  nonlinear  dynamics  of  the  system  and  elasticity  effects  in 
the  towline. 

The  close  proximity  ship  towing  is  under  current  study  by  Office  of  Naval 
Research.  Traditionally,  ship  towing  for  the  military  ships  and  the  others  are  done  with 
long  towlines.  The  Office  of  Naval  Research  is  studying  the  close  proximity  towing  for 
its  innovative  Hull  Programs,  such  as  small  waterplane  area  twin  hull  (SWATH)  ships 
and  their  variations  (such  as  the  SLICE  hull).  During  the  SEA  LANCE  project,  2000,  at 
the  Naval  Postgraduate  School  (Total  Ship  Systems  Engineering  Program)  it  was  shown 
that  one  of  the  main  benefits  of  close  proximity  towing  in  military  applications  is  the 
ability  to  separate  a  combatant  ship  from  a  main  part  of  its  payload. 

In  this  thesis  the  focus  of  our  efforts  is  on  the  nonlinear  dynamics  during  the  loss 
of  directional  stability  in  close  proximity  ship  towing.  Previous  studies  on  the  directional 
stability  of  ship  towing  were  performed  in  1964  by  Abkowitz  who  developed  the 
characteristic  equation  for  single  body  towing,  and  by  Bernitsas  and  others  who 
developed  the  criteria  for  stability  of  Abkowitz  4,h  order  characteristic  equation.  The 
underlying  assumption  in  those  studies  was;  the  leading  ship  is  on  a  constant  straight-line 
motion  of  a  point  mass.  This  is  a  valid  approximation  for  long  towlines  while  it  may  not 
capture  the  combined  dynamics  during  towing  in  close  proximity.  Therefore,  it  is 
possible  that  the  existing  criteria  of  directional  stability  may  be  inadequate.  The 
directional  stability  for  the  coupled  system  where  both  the  maneuvering  dynamics  of 
leading  and  trailing  ships  are  considered  was  studied  by  GOKCE  in  2002.  As  a  result,  the 
complete  system  has  now  a  dimensionality  of  8  and  the  characteristic  order  is  8th  order 
instead  of  4th .  With  the  inclusion  of  the  dynamics  for  the  leading  ship,  it  is  necessary  to 
incorporate  path  control  into  the  system.  For  the  demonstration  purposes,  we  include  a 
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rudder  control  law  for  the  leading  ship  only.  This  is  done  in  terms  of  full  state  feedback  in 
terms  of  the  sway  and  yaw  velocities,  heading  and  lateral  deviation  of  the  leading  ship 
from  its  desired  track.  Control  law  responsiveness  is  quantified  via  its  time  constant, 
small  values  of  the  time  constant  designate  a  responsive  controller  while  larger  values 
designate  a  more  sluggish  design. 

In  this  study  as  a  leading  ship  we  will  use  the  SLICE  vessel,  which  is  105  ft,  and 
1 80  tons,  and  as  a  trailing  ship  the  SWATH  ship  Kaimalino  which  is  89  ft,  47  ft  beam, 
10-12  knots  speed  and  217  tons.  Typical  speeds  for  the  combined  vessels  in  the  tow 
configuration  are  up  to  approximately  15  knots.  The  Slice  is  a  high-speed  variant  of  the 
Swath  technology  with  4  underwater  hulls  instead  of  two.  Attached  to  each  hull  is  a  strut 
that  extends  up  to  support  the  main  body.  Figure  1.  shows  the  profile  view  of  the  Slice. [2] 


The  SSP  (Semi  Submersible  Stable  Platform)  Kaimalino  was  the  world’s  first 
high  performance  open  ocean  Swath  ship.  It  consists  of  two  parallel  torpedo-like  hulls. 
Attached  to  the  hulls  are  two  streamlined  struts.  The  struts  extend  above  the  water  surface 
and  support  the  main  body.  The  Kaimalino  also  has  stabilizing  fins  attached  near  the  aft 
end  of  each  hull.  Figure  2  shows  a  profile  of  the  SSP  Kaimalino.  [2] 
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Figure  2.  The  Profile  of  the  SSP  KAIMALINO  [From:  [2]] 

The  general  approach  adopted  in  this  thesis  is  as  follows: 

1 .  First  we  will  use  the  maneuvering  equations  of  motion  in  the  horizontal  plane  for 
both  leading  and  trailing  ship. 

2.  Coupling  between  the  two  ships  will  be  provided  through  the  towline. 
Hydrodynamic  coupling  arising  from  radiation  and  diffraction  effects  will  be 
neglected.  This  can  be  easily  incorporated  into  the  analysis  once  the  effect  of  such 
hydrodynamic  coupling  on  the  hydrodynamic  coefficients  is  established. 

3.  The  system  will  be  linearized  in  order  to  study  its  stability  properties.  At  the  end 
of  this  part  of  study  identification  of  regions  of  stability  is  expected. 

4.  A  nonlinear  analysis  will  be  perfonned  for  identifying  mechanism  of  loss  of 
stability.  The  response  of  the  system  will  be  analyzed  with  the  help  of  both  the 
linear  and  nonlinear  studies. 

The  hydrodynamic  coefficients  of  both  ships  used  in  this  study  are  provided  from 
reference  [2]. 


3 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


4 


II.  PROBLEM  FORMULATION  AND  STABILITY  ANALYSIS 


A.  EQUATIONS  OF  MOTION 

In  previous  studies,  several  different  values  for  tension  and  for  length  of  the 
connection  between  the  two  ships  were  tried  for  the  stability  analysis  of  the  ship  towing. 
[Gokce,  March  2002]  For  every  instant,  a  zero  eigenvalue  was  found.  This  means  the 
system  cannot  reach  to  stable  condition.  This  zero  eigenvalue  was  attributed  to  the  lack  of 
directional  stabilizing  effects  on  the  leading  ship.  When  the  signs  of  the  towline  restoring 
forces  and  the  moments  are  investigated,  it  can  be  seen  that  towing  force  has  a 
destabilizing  effect  on  the  towing  ship.  On  the  other  hand  it  has  a  stabilizing  effect  on  the 
towed  ship.  The  ships  cannot  show  directional  stability  in  the  horizontal  plane,  therefore 
a  rudder  control  law  must  be  applied  to  the  equations  of  motion. 

The  mathematical  model  of  the  surge,  sway  and  yaw  maneuvering  equations  of 
motion  in  dimensionless  form  is  shown  below.  The  subscript  1  refers  to  towing  (leading) 
ship  and  the  subscript  2  refers  to  towed  (trailing)  ship.  In  these  equations  u  is  the  surge 
velocity,  v  is  the  sway  velocity,  R  is  the  resistance  of  the  vessel  moving  through  body  of 
water,  and  T  is  the  tension  in  the  connection  (rope,  cable,  or  some  other  mechanism) 
between  the  two  towing  ships.  The  rudder  angle  is  denoted  by  8 ,  the  heading  angle  is  vp , 
and  y  is  the  geometric  angle  depicted  in  Figure  3. The  equations  of  motion  are: 


02  -  Tm)v 2  -  Ynr2  =  ~m2r2  +  Yv2v2  +  Yr2r2  -  T  sin {y/2  +  y)  ( 1) 

{Jzi  ~  Y;-2)f2  ~  Ni2v2  —  Nv2y2  +  Nr2r2  —  Txp 2  sin(^2  +  y)  (2) 

K  -  Yn )Vj  -  Yr]rt  =  -/>yi  +  Yv jVj  +  Yrlrt  +  T  sin( if/,  +  y)  +  Ys8  (3) 

(4i  -  Nn  >i  -  =  Avlv,  +  Nrlt\  -  TxpX  sin(^,  +y)  +  N s8  (4) 

where 

5  =  ~kvY\  -  kvv\  -  kA  -  kyy  1  (5) 
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To  study  stability  of  this  system,  the  equations  must  be  linearized.  We  are  dealing 
with  small  heading  angles,  so  the  small  angle  assumption  can  be  made. 

sin  Wi,2=Wi,2  (6) 

cos^12=l  (7) 

Another  assumption  we  can  make  is  for  the  velocities  of  the  ships: 

ul=ii1=\  (8) 

Therefore,  the  linearized  set  of  equations  of  motion 

0 m2  ~  Yvl  )v2  -  Y.:r2  =  -m2r2  +  Yv2v2  +  Yr2r2  -T{y/2  +  y)  (9) 

(/ 72  ~  4-2  —  442  —  4'2V2  4-2*2  ~  p2(l// '2  +  / )  (10) 

K  -  4  )Vi  -  4l*i  =  +  4l*i  +  44,  +  7)  +  YsS  (11) 

(4,  -4iM  -44,  =Nvlv1+Nrlrl -Txpl{ysx  +  y)  +  Ns8  (12) 

The  kinematics  and  the  geometry  relations  of  the  system  are  described  by 
equations  (13)  and  (14). 

y,  =  ux  sin  y/x  +  vx  cos y/x  (13) 

y2  =  u2  sin  y/2  +  v2  cos  y/2  (14) 

If  we  define 

y  =  y 2  Yi  (15) 

then 

y  =  y2-yi  06) 
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These  geometric  relations  are  explained  in  Figure  3. 


When  the  equation  (13)  and  equation  (14)  are  substituted  into  equation  (16),  we 


get 


y  =  u 2  sin  y2  +  v2  cos  y/2  -  ux  sin  -  vx  cos  y/x 
From  the  Figure  3  it  can  be  seen  that; 

v2  +  xp2  sin  y/2  -  (j,  -  xpl  sin  y/x ) 


(17) 


sin  y  - 


l 


(18) 


which  can  be  rewritten  as 


sin  r  =  y  +  y  (. xp2  sin  y/2  +  xpl  sin  y/\)  (19) 

Using  the  small  angle  assumption  once  more  we  get 

y,  =Y\+  U  (2°) 
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y2=y/2+V2 


(21) 


y2+x  2y/2-[y,-x  xy/x) 

sin  Y  = - - - - - 

/ 

The  next  step  is  to  convert  these  equations  into  a  matrix  form. 


^2  —  A2wV2  A2vt-r2  -®2v  (?^2  7 )  (23) 

r2  —  A2rvv2  +  A2n.r2  +  B2r{if/2  +  y)  (24) 

Vl  =  4vvVl  +  Ayr* i  +  5lv  (^1  +  y)  ■ +  CUS  (25) 

A  =  ArvV  1  +  (^1  +  ?)  +  (26) 


The  coefficients  in  equations  (24),  (25),  (25)  and  (26)  are  found  after  the 
following  mathematical  steps: 


4„  =  1(1,2  -  N„  )Y,2  +  (N,2Y„ )] -  l(m2  -  Y„  X/,2  -  JV„ )  -  JV„7„]  (27) 

4.,  =  [(>',,  -MjX/,2  -  AT,2)  +  (Af,.2r,,)]^[(m,  -F„,)(/..2  (28) 

4.  =  H/rf  -  -r„Va  -Nj-N'Jn]  (29) 

A„=[(Y,2Nj+Njm2 -Yi2)]+[(m2-Yi2ilz2 -N„)-Na¥„]  (30) 

4.,  =[fo  -m2)jv«  +4,(m,  -rJM(m2  -r«X42  -  WJ-1M-J  00 

4,  =  [-  TNn  -  Txjm ,  -  7,  )J-[(»7  -  7  , )(/.,  -  AAJA]  (32) 

4„  =[(4, -N„K  +(JVv,n,  )]-[(”>, -rjj„ -i V„)-(V„7„]  (33) 

4. ,  =  [(4  -™,X4i  -^„)+(A',,4)Mk  -n,X4  -4J-V4,,]  (34) 

4,  =  1(4  -  N„  )T  -  (Y„Txp, )]a [(,«,  - ¥„  X4,  -  JV„ ) -  N„Y„  ]  (35) 
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k, =[(n,Jv«)+^.(”>,  -n. )]-[('»,  -n,X/„ -n„)-n„  y,  ,]  w 

Ar  =  [fc  -”■>*,  +  -n, )]-[(»,  -n,X^, -Nn)-NHT„]  (37) 

.  =l77V(,  +3i„(m,  -n,)M(”!,  -4X4  -Af„)-Af„i;,]  (38) 

C„ =[(/-i -Nn)Ys  +  NsYn\^[(mi  -F„X/;1  -JV„)-JV(,yH]  (39) 

C„  =  M„  +  W, (m,  -  Y„ )]  -  [(,»,  -  Y„  )(/;l  -  N„ )  -  N„Y„ ]  (40) 


The  hydrodynamic  coefficients  are  taken  from  reference  [2]  with  the  exception  of 
Ys  and  Ns.  There  is  no  data  available  for  the  rudder  design  of  the  SLICE,  so  we 

assumed  a  ship  turning  radius  of  three  ship  lengths  under  fifteen  degrees  of  rudder,  and 
calculated  the  necessary  values  of  the  rudder  hydrodynamic  coefficients  to  achieve  that 
turning  radius. 

The  control  gains  ky,kv,  kr and  ky  are  calculated  by  standard  pole  placement 
techniques  by  utilizing  the  Matlab  program  shown  in  the  appendices. 

The  system  matrix  of  8x8  for  the  stability  analysis  is  as  follows: 


^2 

V2 

c 

r2 

L 

L 

=  M 

ri 

y2 

y2 

Ti 

y  i 

Yi 

Yi 

where 


(41) 
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M= 


o  o 


b2v 

_B^_ 

B~,  x  o 
B2  +  2v  "2 

^2vXpl 

l 

l 

2v  / 

l 

Blr 

_B2r 

1  B2'XP2 

BlrXpl 

l 

l 

B2'- +  l 

l 

Bw 

_K_ 

Clvky 

B\vXp2 

Blv  +  KXpl 

l 

l 

l 

lv  / 

Blr 

_  B'r  _ 

C  k 

B\rX  pi 

B  +KX>' 

l 

/ 

'-'1  rKy 

l 

K+  l 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

-CiA 


-clrkv 


(42) 


Once  the  matrices  are  formed,  we  can  predict  the  stability  of  the  system  by  an 
eigenvalue  analysis.  This  analysis  is  outlined  in  the  next  section. 


B.  DEGREE  OF  STABILITY 

The  critical  eigenvalue  of  the  system  is  the  one  that  has  largest  real  part.  This 
largest  real  part  is  the  degree  of  the  stability  of  the  system.  If  we  have  a  negative  real  part 
it  means  we  have  stable  system.  On  the  other  hand  the  eigenvalue  with  a  positive  real  part 
indicates  that  the  system  is  unstable. 

C.  REGIONS  OF  STABILITY 

The  results  can  be  summarized  in  a  single  graph  designating  regions  of  stability 
and  instability.  This  graph  is  shown  in  Figure  4,  which  was  produced  by  utilizing  the 
Matlab  code  shown  in  Appendix  A.  The  graph  shows  the  critical  value  of  the  control  time 
constant  for  stability  versus  the  towline  length,  and  is  parameterized  by  the  tension  in  the 
towline.  Combinations  of  (Tc,L,T)  below  the  corresponding  curve  will  yield  stable 
response,  while  combinations  that  are  located  above  the  curve  will  result  in  system 
instability. 
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1.2  1.4  1.6  1.8  2 


Figure  4.  The  Stability  Region 
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III.  NONLINEAR  STABILITY  ANALYSIS 


A.  HOPF  BIFURCATION 

In  all  cases  of  the  loss  of  stability  that  we  defined  in  the  previous  chapter,  a  pair  of 
complex  conjugate  eigenvalues  has  zero  real  parts.  This  is  called  Hopf  Bifurcation  and  is 
accompanied  by  periodic  solutions  (limited  cycles).  The  Hopf  Bifurcation  occurs  when  a 
pair  of  complex  conjugate  eigenvalues  crosses  transversally  the  imaginary  axis.  When 
this  happens,  the  system  will  deviate  from  a  steady-state  solution  in  an  oscillatory 
manner.  This  deviation  can  be  either  supercritical  or  subcritical.  For  the  supercritical 
situation,  stable  limit  cycles  form  after  straight-line  stability  is  lost.  Periodic  solutions 
appear  after  the  initial  loss  of  stability  and  the  resulting  periodic  solutions  are  stable.  For 
the  subcritical  cases,  the  periodic  solutions  appear  before  the  initial  loss  of  stability  and 
the  resulting  periodic  solutions  are  unstable.  In  such  a  case,  the  domain  of  attraction  of 
the  equilibrium  point  becomes  progressively  smaller  as  the  critical  point  is  reached,  and 


c 


(b)  subcritical  bifurcation 


Figure  5.  Supercritical  and  Subcritical  Poincare-Andronov-Hopf 
Bifurcations 
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in  fact  it  shrinks  to  zero  at  the  critical  point.  Therefore,  the  equilibrium  point  is  unstable 
under  finite  disturbances  and  in  practice  it  may  lose  its  stability  even  though  the  linear 
analysis  is  predicting  it  to  be  stable.  Classification  of  the  nature  of  these  bifurcations  is 
the  subject  of  this  chapter.  In  order  to  evaluate  the  stability  of  the  limit  cycles  of  Hopf 
Bifurcations,  a  nonlinear  analysis  is  needed  to  be  performed. 


B.  FORMULATION 

In  the  previous  chapter,  in  order  to  get  the  degrees  of  stability  and  the  regions  of 
stability,  the  equations  of  motions  were  linearized.  In  this  way,  the  linear  stability  of  the 
ship  towing  was  studied.  In  this  chapter  the  nonlinear  dynamics  of  the  system  will  be 


studied. 

The  nonlinear  forms  of  the  equations  were  described  before: 

(m2  ~  yv2)v2  ~  Y-2r2  =  -m2r2  +  Yv2v2  +  Yr2r2  - T sin {y/2  +  y)  (1) 

(7z2  _  ^ rl)^2  ~  Ni2V2  ~  4W2  T  ^ r 2^2  ~  Yx  p2  Sm(^C  T  Y)  (2) 

K  -  4  4  -  4h  =  ~m\rx  +  4vi  +  Yr/\  +  T  sin(  Wi  +  r)+ Yss  (3) 

(4 1  -  NnVi  ~  =  Kivi  +  Nr\r\  ~  TxPi  sin(^i  +  r) + Nss  (4) 

When  we  write  these  equations  in  the  matrix  form;  they  are  as  follows: 

4  =  4vvV2  +  A2vrr 2  +  B2v  sin  {V 2  +  Y)  (43) 

r2  =  4nT2  +  K/i  +  B2r  sin {Vi  +  Y)  (44) 

+  AvA  +  K  sin  +  /)  +  Ch,S  (45) 

K  =  AnA 1  +  A, A  +  K  sin {Yi  +  r)  +  ClrS  (46) 

y  =  u2  sin  y/2  +  v2  cos  y/2-ux  sin  y/x  -  v,  cos  y/x  (47) 
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The  first  step  for  the  nonlinear  dynamics  analysis  is  to  expand  the  y  term  in 


Taylor  series  up  to  the  third  order.  We  have: 

y2+x2sml//2-(yl-x1sml/fl) 

sin  7  = - - - - - - -  (18) 

by  utilizing  thesinv|/  taylor  expansion: 

1  , 

smy/  =  Y —  y  (48) 

6 


y= Ml 


IT+xJ  Vi-J,  Yy  Wx-yY I 


(49) 


This  alpha  value  is  applied  to  the  equations  (43),  (44),  (45),  (46).  In  order  to 
calculate  the  nonlinear  dynamics  coefficient,  we  have  to  form  a  matrix  of  g  values.  By 
these  we  refer  to  the  nonlinear  terms  in  the  equations  of  motion  up  to  third  order.  These 
g  (x)  values  are  obtained  from  the  equations  (50),  (51),  (52),  (53),  (54). 


=  4wV2  +  Av/ 2  +  B2v  ^MYi+T)  (50) 

g\(x) 

A  =  A2,rV2  +  A2r/2  +  +  (5  1 ) 

g2(X ) 

Vi  =  4wh  +  Aurf  i  +  K^{Y\  +  r)  +  ClvS  (52) 

A  =  Ainvi  +  +  K  si n{Yx  +  r)  +  Ch.S  (53) 

g4(x) 
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(54) 


y  =  sin  y/2  +  v2  cos  y/2  -  sin  yx  -  vt  cos  Y\ 

v 

g5(x) 

With  these  g  (x)  values  a  matrix  of  lx7  is  formed: 


g  0)  = 


fiW 

g2(x) 

g3(x) 

g4(x) 

g5(x) 

g6(x) 

g7(x) 


(55) 


Next  step  is  to  apply  the  Taylor  series  expansion  for  the  sin(\|/  +y)  part  of 

equations  (50),  (51),  (52),  (53)  and  keep  the  nonlinear  terms  up  to  third  order.  The  terms 
g6(x)  and  g7(x)  are  equal  to  zero.  In  order  to  account  for  rudder  angle  saturation  we 
substituted  the  discontinuous  saturation  function  with  a  smooth  function  based  on  the 
hyperbolic  tangent  of  the  rudder  angle.  Since  we  are  dealing  with  nonlinear  analysis,  the 
nonlinear  terms  up  to  third  order  are  kept  and  the  following  coefficients  are  obtained. 


1  3  ,.,3  1  ,  1  ..  ..2  1  ..2 


gl  Blv^6L3  Xp'1^1  2 L 


1  2  1  2  1  1  2 

2  y2Wi  +  ^  r3  T1T2  ttttTi  y2  +  —xplw2  „ ,  T2  Ti  Yi 


2  v 


2L 


L 


2L 


-^X,Xp2¥W2  -Jj3Xptxl2ll/y;+—X2p\Xpltfvl  +^yly2Xp2^2  --^y^2Xp2Vl 

1  1  3  1  1  3 

+  -^-yiyixpiV'i  -  — jyiy2xp2yx2  +  -yy Xx  pXx  p2y  w  2  -—yylxplxp2i/sli/s2 

L  O  L  L  O  L 

~  ^ jj,  y \X P\X p2}y y\xpixp2¥\  w2  +  ^ ^  y2xP\X  p2^\  w \  +  ^ ^  y2xPiXp2ysxyf2 


~^y2xplxp2Y3y2  -jyy2xplxp2V \W2  ~-^xpw\  — y2xp2y\ 
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+  j2yiXp2V2+j2yiy2V2-JjTXl^2-^y2V22+^yiV22-Jj2x2p2V'2 

1  3  2  1  2  1  1,1 

+  77TXP1XP2^^2  -7V/l^  -TT yiXpXXp2¥x¥2  +T7iy2XMV2+^XpW2 

oL  L  L  oL  L 

-6F*x^>2  “  J^y2x2p2V2  +  -^y^piYl  +  v22xp2^23  -—ylxplyr2 

-jjTylxpiVl  +  j^ylxplri  +  jjTyixl2V2-jjjylXp2V2  +  Y^Tytxp2V/2 
~^ylxpxV\ +^yxP^ ~^y2xWx  -^xM  +jxpM ~}ixM 

1  3  „3  .  1  3  .  1  2  ,„3„,2  1  3-j 

6F’’1 


,3,, ,3  1  2  .,,2,, ,3  .  1  3  1 


& 2  =  B2r{—^ylxplx  2Y^i2  +—Tx  lx  2yr1  y/2  +  — T y\  -—^-771^2 

JoL  IzL  oL  1L  1L 


+  ^rry\Y12  ~  -^ttx\xYIY2  +  -yrr  y  2X  vxY\y  2  ~  +  tt  v, v2^2 


+  -^2xpxxp2y\yI  -  -jtx Pxx p2y xvl  -  x-y\  -  Yjyy2^2  +  jyyixPixp2YiY2 

-~M}yXXPXXp2V  1V2  -Jj3x2pXXP2YxY2  ~  J^XpXx2p2YxY2  +  X pX^Wl 

--Z^yiXpXXp2YxY\-^y1XXp2Y2+T^y1XXp2Y\-l^y2xlpxY\-l^y12XP2Y2 


+  jjyyxx2P2Y22  +  jYyxxP2Y2  -jjjy22xPxYx  +  -^jy\xPxY\  +  ^jryl*PiYl 


-^yXpxYxYl  -  T7T •y2xp2^22  -  +  yr yxxPxYxY2  ~  T7xPxVl 

ZL  ZL  OL  L  OL 
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1  1  1  3  1  1  3 

~Txp\V^\  ~rxpiW2  ,  r3  y\ y 2X pW i  -ps  y ix p\x piV w i  +  ~7tj y 2X p\x piV^ \V 2 

L  L  oL  L  oL 


+ 


1 

6 U 


y  2X  p\X  p2^  \  V 2 


1 

36  L 


y  2X  p\x  p2^  i  V 2 


1 

YU 


yf  y  2  + 


1 

YU 


y\y\  ~ 


1 

6Z7 


3 

2 


6^3  -V2  2L2  X p2^ 2  +  y ^y 2X p2^ 2 


1  3  1 

~^~jr  y \y  ix  P2V  2  + 


1 

6L 


xp>^ 


3 

1 


1 


21 


1 2Z 


-yUpiVl  + 


21 


3  Vl^pl^l2 


31 


x  p2 


3 

2 


+  -^XpiVlvl  -jry2xp2V2 


=  +  ^f  ~YUyly/l  ~YUyfy/l  ~UUx>l  +  YUy^ 


1  ,2  1 


1 


+  —  v2^r  +  y2xPixP2Vxvl  -  j2y2xPiVi  -  -^xp2v'v2  +  ~^xP2vW\ 


+  jyy^piV'i  -  y^U^xyI  +  jy y^i^  + 


+  ■ 


1  2  3  1  2  3  1  22^  22^2 


12L 


3  J2Xpl^l  +  ,  ~  r3  y 2X p2^ 2  0  r3 


12Z 


1V  y2xP2V2  ~jjyy2xp1wl  -jjry2xP2^2 


+  7T7T y ixl 2V 2  ~  7T7T y  -  -7TyyxxpxxP2¥,¥l  +  T^ryxlxxP2¥x¥l 


2  L 


2 V  '  6L 


12  L 


+  W77yyixplxp2^¥2  +77Txpixp2^>23  -ttxpixp2^>2  -T7yylyixP2V2 

3oL  ol  L  OL 

~7y  y \y  2X  p2^  2  T  j  2  y  2X  p2V^\V^2  T  r3  3^  2X  p\X  p  2}^  1  V  2\X  p2^  2  ,  3  J7 1  X  p2^2 

L  oL  ol  IzL 


1 


Yjryixpi¥i  +  +  -Yjjy^vl 
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1111  1  ,3 

-fjyiXpiVWi + y\yixP\¥\  +-ty&Pww2  ~~^y^xP \xP2ViVi  —^y^Xp^K 


+  7T);iVP2^i^2-T777y2xi,1x;,2^f2 

L  Jo  L 


3 ’"3  — X2 ,  v/]  -  —  X  ,  i//.3  +  — V.  V-2 

2Z2  p  ^  31  ^  ^  2Z3  J 


6  L 


3  X^'  -JZXpl¥2  ~^y2yi  ~^¥'  -^yy^p^iYl  - 


j_ 

2Z3 


1  3  1 

6* 


J_ 

2Z3 


1  1C 

-^y^XpiXpiYlYi -^xPixP2YxY2}-jr{ky,krrM  +kvkyy\Yl  +  Mvvi^i 


+KKr\Y\  +  kJclytyi  +  2£,MvW1  +  krklrX  +  2KKKviri  +  ^-KKK^yWx 


+2KkrkyrxyxYx  +  klkrvY\  +  +  KKri\ + KktvX  +  KkXy\  +ykvY 


2  2 


,2,.  ..2  .  l2l  2 


1 


_3,„3 


+  Tklvl  +~klr\  +~klyl } 


g4  =  v3 ^  K  -^*>2 

+  ^r  y^p^p^^  -JTyiXpiX-^XpiYXi+^XpiXYi+jryiXpiX 


1  2  2-1 


21- 


2  ^2^1^  +  ^  +  ^7r^^»2 


1 

12Z3 


21 


~  y\x  p^Y  ] 


+  ■ 


121 


+-^ylxp2Y32  -  j^y2x2p2Y2  ’ 


12Z 


2/. 


-v2xPi^i  2^  y 2X p2Y 2  "*"  0  r3  y\Xp2Y2 


2  V 


1 


21 


2  1  2  3  1  2 

T-Vl  XP2^2  +  . . 


12Z3  -— t  vrxpl^  +  vfxpl^3  ^X>2 


-^yxXpXxplYW\  +  ^j3x2pixP2Y2Y32  +^ylXplXp2YxY2  +-^2 xpXxplYlv\ 
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1  2  1  3  1  1  3 

x p\x piV^ i  ¥ 2  ^  jj,  y \y 2x p2^ 2  jj,  y \y 2x piW 2  ^ ^2  y  2x  Pi¥ \¥  2 

1,1  1  1  1  1  , 

-^y2xp2yly2+-[jy]y2xplyl  +^ylxp2¥x¥i  ~^y2xPxxPi¥x¥2 --^y{y^pm 
+  ^yiXPlXP2VxV2-^^y2XPXXP2Vlvl-^2Xpx¥l  -JfXP^+J^y'yi 


' ~j xl \¥\  ~  ~ xp 2 ¥2  _  ~t>’i2>’2  -  ~ ¥\  ~~  ~j }!\x P:¥\¥i  ~  ~j x],\x Pi¥\  ¥ 2 


1  1C 

-^Wp2¥l¥2  -7^XplXl2¥l¥l}-^{KkrWl+Kkyyx¥x  +kXVl¥x 


+Kkrr \2¥i  +  Kk\,yW\  +  ^AMfTi  +  krklr\y2\  +  ikJcXm  +  2k  kkv^i 


+2Kkrkvriyi¥i  +  ^AA  +  klk^y^  +  KKkv  1 + KKvx  A + KkXy\  +i:K¥i 


2,.  ,.2  .  l2i  .2 „ 


3  v  1  3,1  3  v2U 


In  order  to  get  g5 ,  a  Taylor  series  expansion  of  sine  and  cosine  is  used. 

1  3 

sin  y~¥~~¥ 

X  1  2 

cos  ^  =  1  —  —  lff~ 

g5=¥2-¥t-^(¥2- ¥] )  +  ¥~y~ -(¥¥2  - y¥x  ) 

The  next  step  is  to  transform  coordinates  to  nonnal  form.  This  is  obtained  by  the 
using  matrix  of  eigenvectors  [m]  at  the  critical  point  as  the  coordinate  transformation 
matrix.  The  new  set  of  coordinates  is  called  [z],  so  in  terms  of  the  original  set  we  have 
[x]=[m][z].  The  only  variables  of  [z]  that  need  to  be  retained  in  the  analysis  of  this 
problem  are  first  two,  the  remaining  variables  are  expressed  in  tenns  of  first  two 
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variables  through  higher  order  relationships  and  do  not  affect  the  third  order  expansions 
that  we  utilized. 


y/]  =  muzl  +  ml2z2 


vi  =  m2izi  +m22z2 


r\  —  mnzi  +  m32z2 


yi=m4lZl  +«42Z2 


If/ ^  m51Zl+m52Z2 


v2=m6]z\+m62z2 


r2=milZl  +  ™12Z2 


y2=m8lz1+mS2z2 


When  the  above  terms  are  substituted  into  the  equation  of  motion,  we  get: 


1  .  3  33  33\  1/32.2 


Al  r3  (W41  mSl  m5lXp2  m\\Xp\)  2Il  Kl^2  +  W41W51  +  W1 1^51^1 ) 


6L 


12  1  i 

-m5i~—m5iXp2 -  —  mnxPi+^{mAxmu-mAXmix-miXm%xxp2  +  mAXmnxp2 


-mfxmnx2pl  -  mnm\xxpX  -  m24lm5lxp2  -  mnm24lxpl  +  m2,m4Xx2pX  -  mum251x2p2xpl 


~m2lm5lx2plxp2)  -  -jj(m25lm& Xxp2  -  m4lm25lxp2  -  m4lm5lm&l  +  mnm5lmslxpl 
-mx  xm4lm5lxpX  +  m,  xmlxxpXxp2 )  -  (m5>81  -  m41m521  +  mx  tm2txpt )  +  (m4]m5tm,txp2 


+mnm4lmslxpl-mnm5lm8lxplxp2+mnm4lm5lxplxp2) 
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A  2  (m42  nhl  m52Xp2  m\2Xp\)  r^j2^tn‘i2X p2  m\2m52Xp\ )  ^  m52  3^  m52Xp2 


^  m\2Xp\  2^3  (m42m$2  m42m82  m52mS2Xp2  fn42fn52Xp2  ^mi2mS2Xpl  m\2m%2X p\ 

-m;2m52xp2  -mnm242xpl+mf2m42x2pl  -ml2m252xplx2p2-mf2m52x2plxp2)-^(mS2m22xp2 
-m42m22xp2  +  2ml2m52  -  m42m52mS2  +  ml2m52mS2xpl  - ml2m52m42xp]  +  mnm252xpixp2) 
-Jj-(m252m%2 -m42m;2  +ml2m22xpl)+^(m42m52mS2xp2  +ml2m42m82xpl  -mnm52muxp]xp2 


+mnm52m42xplxp2) 


A3  B2v{ (,m4\m42  mnm%2  m5im52Xp2  m\\m\2X p\+ m42m%\  m4im«2  m5imS2Xp2 

+m42m25lx2p2  -  2mflm82x2pl  -  mnm2gixpl  -  ml2m24lxpl  +  m2um42x2pl  -  ml2m25lxplx2p2 


mnm52x2xD2)-^^-m2m52x2  -  m4lmum%2  +  m4lm42m8l  +  m5lm52m8lx 


1  \  52  p\  p2 


V  2 


51  52  ^2  *'41  *'81  fr82  1  "Ml"*' 42"*  81  1  "*51"*52"*’81,'v/>2 


-m4lm5lm52x2p2  +2mllml2m8lx2pl  +mnm81m82xpl  +m4lm42m5lxp2+mllm4lm42xpl -mnmum4Xx2pl 
-m4im5in^2xp2 -m42m5in^xp2 -mnm4]n^xfA -m^m^x^  -mnm42m8l  +mum52m%]xp]xp2 
+ml2m5lm%lxplxp2  +  mnm5lm82xplxp2  -  mnm4lm52xplxp2  -  ml2m4lm5lxplxp2 

-mum42m5lxplxp2  +  mum51m52xplx2p2  +  mnml2m5lx2plxp2)--m2lm52  -  —  (m2tmnxp] 


+m; \mS2  -m42m;l  +  mnm2lxpl)-j2(2m5lm52muxp2  +  m5lmH2xp2-2m4lm5lm52xp2 


-m42m5 jxp2  +  m41m42m5l  -  m41m52m81xp2  -  m41mS2mgl  -  m41m5lm82  -  m42m5lmsl 
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+mnm52mgixpi  +  m12m51mglxpl  +  mum5lmS2xpl  -  mnm4lm52xpl  -  ml2m4lm5lxpl 


1 

-mum42m5lxpl  +  2  mnm5lm52xplxp2  +  ml2m5lxplxp2  +  mnmnm5lxpl)--{m5lm52mn 


-m4lm5lm52  +mnm5lm52xpl)-^2;(m24lm52  +  m24Xm52xp2  +  mzum52xpl) 


L\4  ^2v^2Li  42  mSimS2  m5im52Xp2  m \\m \2X p\  +  m 4\m §2  m42m%\ 


-m;2m8lX2p2  +  m4lm;2x2p2  -  m242m5lxp2  -  mnm242xpl  +  m22m4lx2pl  -  mnm22xplx2p2 


1  3 

~mi2m5lXplXp2)  TT^~Zm5im52X d2  m52m %\X nl  5\m 5im%2X b2  ~  ^41^52^ 


L  2 


-2 m42m5lm52xp2  -  m4lm52m&2  -  m42m52mu  -  m42m51m82  +  ml2m52mslxpl 
+mnm52mS2xpl  +  ml2m5lmS2xpl  -  ml2m4lm52xpl  -  mum42m52xpl  -  ml2m42m5lxpl 


+  mnm22xplxp2  -  2 ml2mslm52xplxp2  +  mnml2m52x2pl  +  m22m5lx2pl)  - -m5lm22 


-j(2m5lm22xp2  +  m5lm52mS2  -  m42m5lm52  +  ml2m5lm52x pl)  -  jj-(mnmf2xpl 

+ml2mu  -  m4lm22  +  mnm252xpl)  +  —  (m42mHlmH2  -  m4lm42mS2  -  m5Xm52m%2x2p2 

+m42m5lm52x2p2  -  m4ltn42m52xp2  -  m12m4lm42xpl  +  mum12m42x2pl  +  m4lm52mS2xp2 
+m42m52mslxp2  +  m42m5lmS2xp2  +  ml2m4lmS2xpl  +  ml2m42mslxpl  +  mnm42mt2xpl 
-mnm52mS2xplxp2  -  mnm52mt2xplxp2  -  ml2m51mS2xplxp2  +  ml2m4lm52xplxp2 


+mum52m42xpXxp2+mX2m42m5XxpXxp2-mX2m5lm52xplx2p2-mnmnm52x2pXxp2) 


23 


Following  the  above  background  calculations  we  are  now  ready  to  summarize  the 
nonlinear  analysis  for  our  system. 


C.  THE  CONTROL  LAW 

A  typical  full  state  feedback  control  law  will  be  utilized  for  the  nonlinear  analysis 
again.  In  the  below  equation  So  is  the  feedback  rudder  angle  for  small  deviations. 

8  =  k  ys  +  Ay  v, +  Ayr, +  k  y,  (56) 

In  order  to  capture  the  effect  of  rudder  saturation,  the  control  law  is  modeled  as: 


8  =  8  sat  tanh 


\0  sat  J 


(57) 


The  hyperbolic  tangent  function  is  used  instead  of  hard  saturation  function 
because  of  its  analyticity  properties. 

Utilizing  the  taylor  expansion  again; 


1 

sin  y/  =  y/-  —  y/ 


3 


8  -  8o~ 


3  8 


Si 


(58) 


In  the  Matlab  (Appendix  A  and  B)  code  for  the  nonlinear  analysis  saturation  limit 
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on  the  8  is  set  to  0.4. 


D.  INTEGRAL  AVERAGING 

When  the  equations  of  motion  are  written  in  the  compact  form, 

X  =  f{x),  X  =  [v2,r2,vl,i\,y2,yl,y/2,y/lf  (59) 

The  system  can  be  written  in  the  form 

X  =  AX  +  g(x)  (60) 


A  is  the  Jacobian  matrix  of  f  (x),  and  g  (x)  contains  all  nonlinear  term  of 
equations  (50),  (51),  (52),  (53),  (54). 

If  we  introduce  the  transformation  matrix  (T)  of  eigenvectors  of  A 


T  1  =  [m;j] 

i,  j  =  1,2, 3, 4, 5, 6, 7, 8 

(61) 

T  *  =  [nij] 

i,  j  =  1,2, 3, 4, 5, 6, 7, 8 

(62) 

linear  change  of  coordinates 

x  =  Tz 

z  =  T  lx 

(63) 

will  transform  system  (X)  into  its  normal  form 

z  =  T-1ATz  +  T~1g(Tz )  (64) 

with  this  transformation  we  get 

O-wO  0  0  000 

-wOO  0  0  000 

0  OCjOOOOO 
0  00  C7  0000 

(65) 

0000  C3  000  v  ' 

0  000  0C400 

0  00  0  00  C50 

0  000  000  c6 
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The  expressions  for  g  are  given  by: 


hlZl  hiz2  hiZlZ2  KaZ\Z2 

^2lZ\  T  ^22Z2  T  ^22Z\  Z 2  ^24Z\Z2 

g  —  l31Z  j  +  /32z2  ^33^1  Z2  ^  24Z  \Z  2 

^  41  Z 1  ^42Z2  1 42  Z 1  Z  2  ^44Z1Z2 

_hlZl  ^  52Z  2  +  ^  53Z 1  Z  2  ^54  Z  \Z  2  J 

(66) 

If  we  substitute  equation  (63)  into  (60),  we  get 

z,  =  -wz2  +  q  jzf  +  r12zf  z2  +  r^z’  +  r14z2 

(67) 

Z2  =  — WZj  +  r2iZj  +  7"22Zj  Z2  +  r23ZjZ2  +  1 24Z 2 

(68) 

These  equations  can  be  written  as, 

z  i  —  —  wz2  +  Tj  (Z| ,  z2 ), 

(69) 

Z2  —  — WZj  +  F2  (zi  j  ^2)’ 

(70) 

Introducing  the  polar  coordinates  for  the  equations  (69)  and  (70), 

Zj  =  R  cos  0,  z2  =  7?  sin  # 

we  get 

(71) 

7?  =  Fj  (7?,  #)  cos  0  +  F2  (7?,  6?)  sin  6 

(72) 

1 10  =  wR  +  F2(R,d)cos6  -  Fl(R,d)smd 

(73) 

Equation  (73)  gives 

R  =  P(0)R3 

(74) 
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The  function  P{0)  is  a  2  n  periodic  function  in  the  angular  coordinate  6 .  In 
order  to  get  the  equation  with  constant  coefficients,  equation  (74)  is  averaged  over  one 
cycle  in  6 .  This  way  we  get 


R  =  KRi 


and 


(75) 


K  = 


—  \Kp{d).dd 
In  •'° 


(76) 


After  making  the  integral  calculations  we  get  the  simplified  form  of  the  equation  (76) 


K  =  ~\3rn+ru+r22+3r24] 


(77) 


The  indicated  coefficient  are  given  by, 


hi  —  nuh\  nn^2\  n\ih\  PJ41  n\sh\ 

hi  ~  niJl4  n\2^24  ni3^24  °14^44  ^15^54 

r22  ~  n2\hi  n22^22  n  23^33  n 24^43  n  25^53 

r24  ~  n2^\2  n22^22  ^23^32  ^24^42  °25^52 


(78) 


The  coefficients  nti  are  the  elements  of  the  inverse  of  transformation  matrix  T, 

which  was  defined  as  the  matrix  of  eigenvectors  of  the  linearized  matrix  A. 

Results  in  tenns  of  the  nonlinear  stability  coefficient  K  are  shown  in  Figure  6.  It 
can  be  seen  that  the  short  towline  lengths  result  in  supercritical  bifurcations,  while  larger 
lengths  result  in  subcritical  bifurcations.  The  implications  of  this  on  system  stability 
under  finite  disturbances  were  discussed  previously.  It  should  be  pointed  out  that  the 
transition  from  supercritical  to  subcritical  bifurcations  as  towline  length  is  increased  is 
dependent  on  the  hydrodynamic  characteristics  of  the  two  ships  and  is  not  necessarily 
true  for  all  ships.  The  procedure  that  we  have  established  in  this  work  should  be  followed 
in  order  to  assess  the  type  of  corresponding  bifurcation. 
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Kcoef 


L 


Figure  6.  Nonlinear  Dynamics  Coefficient  vs.  Towline  Length 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

This  thesis  presented  a  comprehensive  study  of  the  stability  region  properties  and 
the  nonlinear  dynamics  response  of  the  two-body  ship-towing  problem.  The  results  of 
this  study  can  help  in  establishing  rational  guidelines,  which  should  be  followed  in  order 
to  ensure  stability  and  safety  of  close  proximity  ship  towing  operations.  The  main  results 
from  this  study  can  be  summarized  as  follows: 

1 .  The  traditional  towing  stability  model  cannot  predict  the  dynamics  of  towed  ships 
in  close  proximity.  One  of  the  reasons  of  this  result  is  excluding  the  dynamics  of 
the  leading  ship  in  the  overall  analysis. 

2.  Path  stability  depends  on  classical  parameters  such  as  tension  and  towline  length, 
but  also  on  additional  terms  such  as  the  leading  ship’s  control  time  constant. 

3.  For  very  large  values  of  the  towline  length,  the  effect  of  the  leading  ship  control 
becomes  less  pronounced.  Therefore,  previously  reported  results  in  the  literature 
can  be  considered  as  a  large  towline  length  of  the  two-body  model  developed  in 
this  work. 

4.  The  primary  mechanism  of  loss  of  stability  is  in  terms  of  Hopf  Bifurcations  to 
periodic  solutions.  The  bifurcation  can  be  classified  as  either  supercritical  or 
subcritical  in  terms  of  a  compatible  nonlinear  stability  coefficient. 

5.  Hopf  Bifurcations  tend  to  be  supercritical  for  small  values  of  towline  length, 
while  they  tend  to  turn  into  subcritical  as  the  towline  length  is  increased. 

B.  RECOMMENDATIONS 

Recommendations  for  continuing  studies  are  the  following: 

1 .  Analyze  the  properties  of  stability  and  bifurcations  including  the  towline  elasticity 
into  the  model. 
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2.  Analyze  the  stability  properties  of  different  vessels  and  investigate  the  occurrence 
of  static  bifurcations. 
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APPENDIX  A.  MATLAB  CODE  TO  FIND  THE  REGION  OF  STABILITY 


%  Stability  graph  -  TC  vs .  L  -  constant  T 

"o 

%  Parameters: 

%  T  =  Nondimensional  tension 

%  L  =  Towline  length 

%  TC  =  Control  law  time  constant 

%  xpl  =  Attachment  point  (leading  ship  -  positive  forward 

of  amidships) 

%  xp2  =  Attachment  point  (trailing  ship  -  positive  aft  of 

amidships ) 

"6 

%  Constants 


ul 

= 

i; 

u2 

= 

i; 

ml 

= 

0 .018078; 

m2 

= 

0.018; 

Izl 

= 

0 . 0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0 . 07893; 

Yv2 

= 

-0 .1183; 

Yrl 

= 

-0 . 004044; 

Yr2 

= 

-0 . 0042; 

Nvl 

= 

-0 .016428; 

Nv2 

= 

-0 . 0187; 

Nr  1 

= 

-0 . 010332; 

Nr2 

= 

-0 . 0176; 

Yvdotl 

= 

-0 . 051328; 

Yvdot2 

= 

-0.0184; 

Yrdotl 

= 

0 . 005617; 

Yrdot2 

= 

-0 . 0011; 

Nvdotl 

= 

-0 . 001945; 

Nvdot2 

= 

-0 . 0008489; 

Nrdotl 

= 

-0 . 00564; 

Nrdot2 

= 

-0 . 0090; 

Ydel 

= 

0 .0103; 

Ndel 

= 

-0 . 0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0 .5; 

xp2 

= 

0 .5; 

"o 

%  Enter  constant  T 


T  =  input ( ' T  =  ' ) ; 

"5 

%  Start  loop  on  length 

'o 

for  iL  =0.1:0.01:2.0; 

indexL  =  indexL+1; 

indexTC  =  0; 

L  =  iL 

L_v (indexL) =  L; 

A3  =  1/L; 


31 


B3  =  xpl/L; 

C3  =  xp2/L; 

D3  =  - 1 / L ; 

o, 

o 

%  Loop  on  TC 

o, 

o 

for  iTC  =0.1:0.01:2.0; 
indexTC=indexTC+l ; 

TC  =iTC ; 

bpole  =  -1/TC; 

TC_v ( indexTC ) =TC; 

o, 

o 


%  Setup  the  matrix  coefficients 

o, 

o 


A2vv 

=  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ] ; 

A2vr  = 

[ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) )  +  (Nr2*Yrdot2)  ' 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ]  ; 

B2v 

=  [- ( (Iz2-Nrdot2) *T) - (Yrdot2*T*xp2 ) ^ 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ]  ; 

A2rv 

=  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ] ; 

A2rr 

:  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ' 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ] ; 

B2r 

=  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ^ 

]/ 

[ ( (m2- 

Yvdot2 ) 

* (Iz2-Nrdot2) ) 

- (Nvdot2*Yrdot2 ) ] ; 

Alvv 

=  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl*Yrdotl ) ] ; 

Alvr 

=  [  ( Yrl-ml*ul ) * (Izl-Nrdotl )  +  (Nrl*Yrdot 1 ) ’ 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl *Yrdotl ) ] ; 

Blv 

[ ( (Izl-Nrdotl) *T ) - (Yrdotl*T*xpl) ; 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl*Yrdotl ) ] ; 

Civ 

=  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ' 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl*Yrdotl ) ] ; 

Alrv 

=  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl*Yrdotl ) ] ; 

Alrr 

:  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ’ 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl*Yrdotl ) ] ; 

Blr 

=  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

- (Nvdotl *Yrdotl ) ] ; 

Clr 

=  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ’ 

]/ 

[ ( (ml- 

Yvdotl ) 

* (Izl-Nrdotl) ) 

o. 

- (Nvdotl *Yrdotl ) ] ; 

o 

%  Find 

o. 

control  gains 

o 

c 

=  zeros  ( 4 ) ; 

C  ( 1 , 3 ) 

=  i; 

C (2, 2) 

=  Alvv; 

C (2, 3) 

=  Alvr; 

C  (3,  2) 

=  Alrv; 

C  (3,  3) 

=  Alrr; 

C  ( 4 , 1 ) 

=  i; 

C  ( 4 , 2  ) 

=  1; 

D 

=  zeros  (4,1); 

D (2, 1) 

=  Civ; 

D  (3, 1) 

=  Clr; 
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poles  =  [bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 
k  =  place (C, D, poles ) ; 

Kp  s  i  =  k  ( 1 ,  1 )  ; 

Kv  =  k  ( 1 , 2  )  ; 

Kr  =  k  (1, 3)  ; 

Ky  =  k  (1, 4)  ; 

o, 

o 

%  A  matrix 

o, 

o 

A  =  zeros ( 8 ) ; 

A  ( 1 , 1 )  =  A2vv; 

A  ( 1 , 2 )  =  A2vr ; 

A (1, 5)  =  B2v/L; 

A  (1, 6)  =  -B2v/L; 

A(l,7)  =  B2v+ ( (B2v*xp2 ) /L) ; 

A(l,8)  =  (B2v*xpl)/L; 

A  (2, 1)  =  A2rv; 

A  (2, 2)  =  A2rr; 

A  (2, 5)  =  B2r/L; 

A  (2, 6)  =  -B2r/L; 

A(2,7)  =  B2r+ ( (B2r*xp2 ) /L)  ; 

A  (2, 8)  =  (B2r*xpl)/L; 

A  (3, 3)  =  Alvv- (Clv*Kv) ; 

A  (3, 4)  =  Alvr- (Clv*Kr)  ; 

A (3, 5)  =  Blv/L; 

A  (3, 6)  =  -Blv/L- (Clv*Ky) ; 

A(3,7)  =  (Blv*xp2)/L; 

A(3,8)  =  Blv+ ( (Blv*xpl ) /L) - (Clv*Kpsi ) ; 

A  ( 4 , 3 )  =  Alrv- (Clr*Kv) ; 

A  ( 4 , 4 )  =  Alrr- (Clr*Kr )  ; 

A  ( 4 , 5 )  =  Blr/L; 

A  ( 4 , 6 )  =  -Blr/L- (Clr*Ky)  ; 

A  ( 4 , 7 )  =  (Blr*xp2)/L; 

A  ( 4 , 8 )  =  Blr+ ( (Blr*xpl) /L) - (Clr*Kpsi)  ; 

A  (5,  1 )  =  1; 

A (5, 7)  =  1; 

A  (6,  3)  =  1; 

A  (  6,  8  )  =  1; 

A  ( 7 , 2  )  =  1; 

A  ( 8 , 4  )  =  1; 

"6 

%  Find  eigenvalues 

o, 

o 

B=eig (A) ; 

F=real (B) ; 

H_new=max (F) ; 

o, 

o 

%  Detect  if  H  changes  its  sign 

o, 

o 

if  indexTC==l 

H_old=H_new; 

else 

PROD=H_new*H_old; 
if  PROD  >  0 

o, 

o 

%  H  does  not  change  sign  -  keep  going 
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H_old=H_new; 

else 

o, 

o 

%  H  changed  its  sign  -  find  critical  TC  by  linear 

interpolation 

o, 

o 

index=index+l ; 

TC_crit (index) =- (H_old*TC_v (indexTC) - 
H_new*TC_v (indexTC-1 ) ) / (H_new-H_old) ; 

H_old=H_new; 

end 

end 

end 

end 

"6 

%  Plot  (L,  crit  TC)  curve 
"6 

%plot (L_v, TC_crit ) , xlabel ( ' L ' ) , y label ( ' { T_C }_{ \rm 
crit } ' ) , title ( ' Constant  T ' ) , grid 
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APPENDIX  B.  MATLAB  CODE  TO  FIND  THE  NONLINEAR 
DYNAMICS  COEFFICIENT 


%  Nonlinear  Dynamics  coefficient  graph  -  K  vs .  L  -  constant  T 


%  Parameters: 
%  T 
%  L 
%  TC 
%  xpl 
amidships ) 

%  xp2  = 

amidships ) 


Nondimensional  tension 

Towline  length 

Control  law  time  constant 

Attachment  point  (leading  ship  -  positive  forward  of 
Attachment  point  (trailing  ship  -  positive  aft  of 


%  Constants 


ul 

= 

i; 

u2 

= 

l; 

ml 

= 

0 . 018078; 

m2 

= 

0.018; 

Izl 

= 

0 . 0007; 

Iz2 

= 

0 . 00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0 .1183; 

Yrl 

= 

-0 . 004044; 

Yr2 

= 

-0 . 0042; 

Nvl 

= 

-0 .016428; 

Nv2 

= 

-0.0187; 

Nr  1 

= 

-0 .010332; 

Nr2 

= 

-0 .0176; 

Yvdotl 

= 

-0 . 051328; 

Yvdot2 

= 

-0 .0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0 . 0011; 

Nvdotl 

= 

-0 . 001945; 

Nvdot2 

= 

-0 . 0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0.0090; 

Ydel 

= 

0 .0103; 

Ndel 

= 

-0 . 0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0 .5; 

xp2 

= 

0 .5; 

deltasat: 

o 

ll 

"6 

%  Enter  i 

constant  T 

o 

T 

= 

input ( ' T  = 

%  Start  loop  on  length 


for  iL  =1.53:0.01:2.0; 
indexL  =  indexL+1; 
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indexTC  = 

0; 

L 

iL 

L_v ( indexL) = 

L; 

A3 

1/L; 

B3 

xpl/L; 

C3 

xp2/L; 

D3 

-1/L; 

%  Loop  on  TC 

for  iTC  =0.1:0.005:2.0; 
indexTC=indexTC+l ; 

TC  =iTC ; 

bpole  =  -1/TC; 

TC_v ( indexTC ) =TC; 

%  Setup  the  matrix  coefficients 

o, 

o 

A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2v  =  [- ( ( Iz2-Nrdot2 ) *T) - (Yrdot2*T*xp2 ) ] / [ ( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) - 
(Nvdot2*Yrdot2 ) ]  ; 

A2r v  =  [ (Yv2*Nvdot2 ) + (Nv2* (m2-Yvdot2 ))]/[( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Blv  =  [ ( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

%  Find  control  gains 

o, 

o 


C 

=  zeros ( 4 )  ; 

C  ( 1 , 3 ) 

=  i; 

C (2, 2) 

=  Alvv; 

C (2, 3) 

=  Alvr; 

C (3, 2) 

=  Alrv; 

C (3, 3) 

=  Alrr; 

C  ( 4 , 1 ) 

=  i; 

C  ( 4 , 2  ) 

=  1; 

D 

=  zeros (4,1); 

D (2, 1) 

=  Civ; 
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D  (3,  1) 

= 

Clr; 

poles 

= 

[bpole  bpole-0.05 

k 

= 

place (C, D,  poles )  ; 

kpsi 

= 

k  ( 1 , 1 )  ; 

kv 

= 

k  ( 1 , 2  )  ; 

kr 

= 

k  ( 1 , 3 )  ; 

ky 

= 

k  ( 1 , 4  )  ; 

%  A  matrix 

o 

A 

= 

zeros ( 8 ) ; 

A  ( 1 , 1 ) 

= 

A2vv; 

A  ( 1 , 2  ) 

= 

A2vr ; 

A  ( 1 , 5 ) 

= 

B2v/L; 

A  ( 1 ,  6) 

= 

-B2v/L; 

A  ( 1 , 7  ) 

= 

B2v+ ( (B2v*xp2) /L) 

A  ( 1 , 8  ) 

= 

(B2v*xpl ) /L; 

A  (2 , 1 ) 

= 

A2rv; 

A  (2 , 2  ) 

= 

A2rr ; 

A (2, 5) 

= 

B2r/L; 

A (2, 6) 

= 

-B2r/L; 

A  (2 , 7  ) 

= 

B2r+ ( (B2r*xp2) /L) 

A  (2 , 8  ) 

= 

(B2r*xpl ) /L; 

A  (3,  3) 

= 

Alvv- (Clv*kv) ; 

A (3, 4 ) 

= 

Alvr- (Clv*kr) ; 

A  (3,  5) 

= 

Blv/L; 

A  (3,  6) 

= 

-Blv/L- (Clv*ky) ; 

A (3, 7 ) 

= 

(Blv*xp2 ) /L; 

A (3, 8 ) 

= 

Blv+ ( (Blv*xpl) /L) 

A  ( 4 , 3 ) 

= 

Alrv- (Clr*kv) ; 

A  ( 4 , 4  ) 

= 

Alrr-  (Clr*kr ) ; 

A  ( 4 , 5 ) 

= 

Blr/L; 

A  ( 4 ,  6) 

= 

-Blr/L- (Clr*ky ) ; 

A  ( 4 , 7  ) 

= 

(Blr*xp2 ) /L; 

A  ( 4 , 8  ) 

= 

Blr+ ( (Blr*xpl) /L) 

A (5, 1 ) 

= 

i; 

A (5, 7 ) 

= 

i; 

A  (  6,  3 ) 

= 

i; 

A  (  6,  8  ) 

= 

i; 

A  ( 7 , 2  ) 

= 

i; 

A  ( 8 , 4  ) 

= 

i; 

%  Find  eigenvalues 
"6 

B=eig (A) ; 

F=real (B) ; 

H_new=max (F) ; 

'o 

%  Detect  if  H  changes  its  sign 


if  indexTC==l 

H_old=H_new; 

else 

PROD=H_new*H_old; 
if  PROD  >  0 
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%  H  does  not  change  sign  -  keep  going 

"o 

H_old=H_new; 

else 

"6 

%  H  changed  its  sign  -  find  critical  TC  by  linear  interpolation 
"6 

index=index+l ; 

TC_crit (index) =- (H_old*TC_v (indexTC) -H_new*TC_v (indexTC-1) ) / (H_new- 
H_old)  ; 

H_old=H_new; 

"6 

bpole_crit  =  -l/TC_crit (index)  ; 

poles  =  [bpole_crit  bpole_crit-0 . 05  bpole_crit-0 . 10  bpole_crit-0 . 15 ]  ; 

k  =  place (C, D, poles ) ; 

kpsi  =  k ( 1 , 1 ) ; 

kv  =  k  ( 1 , 2  )  ; 

kr  =  k (1, 3) ; 

ky  =  k  (1, 4)  ; 

A_crit  =  zeros (8); 

A_crit(l,l)  =  A2vv; 

A_crit(l,2)  =  A2vr; 

A_crit(l,5)  =  B2v/L; 

A_crit(l,6)  =  -B2v/L; 

A_crit(l,7)  =  B2v+ ( (B2v*xp2 ) /L) ; 

A_crit(l,8)  =  (B2v*xpl)/L; 

A_crit(2,l)  =  A2rv; 

A_crit(2,2)  =  A2rr; 

A_crit(2,5)  =  B2r/L; 

A_crit(2,6)  =  -B2r/L; 

A_crit(2,7)  =  B2r+ ( (B2r*xp2 ) /L) ; 

A_crit(2,8)  =  (B2r*xpl)/L; 

A_crit(3,3)  =  Alvv- (Clv*kv)  ; 

A_crit(3,4)  =  Alvr- (Clv*kr ) ; 

A_crit(3,5)  =  Blv/L; 

A_crit(3,6)  =  -Blv/L- (Clv*ky) ; 

A_crit(3,7)  =  (Blv*xp2)/L; 

A_crit(3,8)  =  Blv+ ( (Blv*xpl ) /L) - (Clv*kpsi ) ; 

A_crit(4,3)  =  Alrv- (Clr*kv) ; 

A_crit(4,4)  =  Alrr- (Clr*kr ) ; 

A_crit(4,5)  =  Blr/L; 

A_crit(4,6)  =  -Blr/L- (Clr*ky) ; 

A_crit(4,7)  =  (Blr*xp2)/L; 

A_crit(4,8)  =  Blr+ ( (Blr*xpl ) /L) - (Clr*kpsi ) ; 

A_crit (5,1)  =  1; 

A_crit (5,7)  =  1; 

A_crit (6,3)  =  1; 

A_crit (6,8)  =  1 ; 

A_crit (7,2)  =  1 ; 

A_crit (8,4)  =  1 ; 

"6 

%Repeated  coloumns  are  4-5  and  6-7  of  the  Devcs 
J=eig (A_crit) ; 

[D_evcs, V_evls] =eig (A_crit) ; 
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T_matrixl ( : , 1 ) =D_evcs ( : , 1 ) ; 

T_matrixl ( : , 2 ) =D_evcs ( :  ,  2 ) ; 

T_matrixl ( : , 3) =D_evcs ( :  ,  3) ; 

T_matrixl ( : , 4 ) =real (D_evcs ( : , 4 ) ) ; 

T_matrixl ( : , 5 ) =imag (D_evcs ( : , 4 ) ) ; 

T_matrixl ( : , 6) =real (D_evcs ( : , 6) ) ; 

T_matrixl ( : , 7 ) =imag (D_evcs ( :  ,  6 ) ) ; 

T_matrixl ( : , 8 ) =D_evcs ( : , 8 ) ; 

"6 

T_matrix ( : , 1 ) =T_matrixl ( : , 6 ) ; 

T_matrix ( : , 2 ) =T_matrixl ( : , 7 ) ; 

T_matrix ( : , 3 ) =T_matrixl ( : , 1 ) ; 

T_matrix ( : , 4 ) =T_matrixl ( :  ,  2 ) ; 

T_matrix ( : , 5) =T_matrixl ( : , 3) ; 

T_matrix ( : , 6 ) =T_matrixl ( : , 4 ) ; 

T_matrix ( : , 7 ) =T_matrixl ( : , 5 ) ; 

T_matrix ( : , 8 ) =T_matrixl ( : , 8 ) ; 

"6 

ml l=T_matrix (1,1) ; ml2=T_matrix (1,2) ; 
m21=T_matrix (2,1) ; m22=T_matrix (2,2) ; 
m31=T_matrix (3,1) ; m32=T_matrix (3,2)  ; 
m4 l=T_matrix (4,1) ; m42=T_matrix (4,2) ; 
m51=T_matrix (5,1) ; m52=T_matrix (5,2) ; 
m81=T_matrix (8,1) ; m82=T_matrix (8,2) ; 
m61=T_matrix (6,1) ; m62=T_matrix (6,2) ; 

"6 

T_inverse=inv (T_matrix) ; 

N1 l=T_inverse (1,1) ; 

N12=T_inverse (1,2)  ; 

N13=T_inverse (1,3)  ; 

N14=T_inverse  (1, 4) ; 

N15=T_inverse  (1,5) ; 

N1 l=T_inverse  (1, 1) ; 

N21=T_inverse  (2, 1) ; 

N22=T_inverse  (2,2) ; 

N23=T_inverse (2,3) ; 

N24=T_inverse  (2, 4) ; 

N25=T_inverse  (2,5) ; 

"6 

M=inv (T_matrix) *A_crit*T_matrix; 

LI 1_1  =  l/6*B2v/LA3*m41A3-l/6*B2v/LA3*m81A3-l/ (6*LA3) *B2v*m51A3*xp2A3; 
Lll_2  =  -1/ (2*LA2 ) *B2v*m51A3*xp2A2-l/6*B2v*m51A3-2/ (3*L) *B2v*m51 A3*xp2; 
Lll_3  =  -1/ (6*L) *B2v*mllA3*xpl- 

/ (6*LA3) *B2v*mll A3*xpl A3+l/ (2*LA3) *B2v*m41*m81A2; 

Lll_4  =-1/ (2*LA3) *B2v*m41 A2*m81-l/LA2*B2v*m81*m51 A2*xp2- 
1/ (2*L) *B2v*m81*m51A2; 

LI 1_5=  1/ (2*LA3) *B2v*m81*m51 A2*xp2 A2+l/LA2*B2v*m41*m51 A2*xp2 
+1/ (2*L) *B2v*m41*m51A2; 

Lll_6  =+1/ (2*LA3) *B2v*m41*m51A2*xp2 A2-l/ (2*LA3) *B2v*xpl A2*m81*mll A2; 
Lll_7  =  -1/ (2*LA3) *B2v*xpl A2*m81*mll A2-l/ (2*LA3) *B2v*xpl*m81 A2*mll; 

LI 1_8  =  -1/ (2*LA2) *B2v*m41A2*m51-l/ (2*LA3) *B2v*m41A2*m51*xp2; 

Lll_9  =-1/ (2*LA3) *B2v*xpl*m41 A2*mll+1/ (2*LA3) *B2v*xplA2*m41*mllA2; 
Lll_10  =+l/LA3*B2v*xp2*m41*m81*m51+l/LA2*B2v*m41*m81*m51 
+l/LA3*B2v*m41*m81*mll*xpl ; 
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LI 1_1 1  =-l/LA3*B2v*m81*mll*m51*xp2*xpl-l/LA2*B2v*m81*mll*m51*xpl ; 

Lll_12  =+l/LA3*B2v*m41*mll*m51*xp2*xpl+l/LA2*B2v*m41*mll*m51*xpl ; 

Lll_13  =-1/ (2*LA3) *B2v*mll*m51A2*xp2A2*xpl-l/ (2*L) *B2v*mll*m51A2*xpl; 

LI 1_1 4  =-l/LA2*B2v*mll*m51 A2*xp2*xpl-1/ (2*LA3) *B2v*mllA2*m51*xp2*xplA2; 
Lll_15  =-1/ (2*LA2 ) *B2v*mllA2*m51*xplA2 ; 

LI 1  =  L11_1+L11_2+L11_3+L11_4+L11_5+L11_6+L11_7+L11_8+L11_9+L11_10 . . . 
+L11_11+L11_12+L11_13+L11_14+L11_15; 

L12_l  =  l/6*B2v/LA3*m42A3-l/6*B2v/LA3*m82A3-l/ (6*LA3) *B2v*m52A3*xp2A3; 

L12_2  =  -1/ (2*LA2 ) *B2v*m52A3*xp2A2-l/6*B2v*m52A3-2/ (3*L) *B2v*m52 A3*xp2; 
L12_3  =  -1/ (6*L) *B2v*ml2A3*xpl- 

1/ (6*LA3) *B2v*ml2 A3*xpl A3+l/ (2*LA3) *B2v*m42*m82A2; 

L12_4  =  -1/ (2*LA3) *B2v*m42A2*m82-l/LA2*B2v*m82*m52A2*xp2- 
1/ (2*L) *B2v*m82*m52 A2 ; 

L12_5  =  -1/ (2*LA3) *B2v*m82*m52A2*xp2A2+l/LA2*B2v*m42*m52A2*xp2; 

L12_6  =  +1/ (2*L) *B2v*m42*m52A2+l/ (2*LA3) *B2v*m42*m52 A2*xp2 A2 

L12_7  =  -1/ (2*LA3) *B2v*xplA2*m82*ml2A2-l/ (2*LA3) *B2v*xpl A2*m82*ml2 A2; 

L12_8  =  -1/ (2*LA3) *B2v*xpl*m82A2*ml2-l/ (2*LA2) *B2v*m42 A2*m52; 

L12_9  =  -1/ (2*LA3) *B2v*m42 A2*m52*xp2-1/ (2*LA3) *B2v*xpl*m42 A2*ml2; 

L12_10  =  +1/ (2*LA3) *B2v*xpl A2*m42*ml2 A2+l/LA3*B2v*xp2*m42*m82*m52 ; 

L12_l 1  =  +l/LA2*B2v*m42*m82*m52+l/LA3*B2v*m42*m82*ml2*xpl ; 

L12_12  =  -l/LA3*B2v*m82*ml2*m52*xp2*xpl-l/LA2*B2v*m82*ml2*m52*xpl ; 

L12_13  =  +l/LA3*B2v*m42*ml2*m52*xp2*xpl+l/LA2*B2v*m42*ml2*m52*xpl ; 

L12_14  =  -1/ (2*LA3) *B2v*ml2*m52 A2*xp2 A2*xpl-1/ (2*L) *B2v*ml2*m52 A2*xpl; 
L12_15  =  -l/LA2*B2v*ml2*m52A2*xp2*xpl-l/ (2*LA3) *B2v*ml2 A2*m52*xp2*xpl A2; 
L12_l 6  =  -1/ (2*LA2) *B2v*ml2A2*m52*xplA2; 

L12  =  L12_l+L12_2+L12_3+L12_4+L12_5+L12_6+L12_7+L12_8+L12_9+L12_10 . . . 

+L12_11+L12_12+L12_13+L12_14+L12_15+L12_16; 

L13_l=l/2*B2v/LA3*m41A2*m42-l/2*B2v/LA3*m81A2*m82 
-1/ (2*LA3) *B2v*m51 A2*m52*xp2 A3 ; 

L13_2=-3/ (2*LA3) *B2v*m51A2*m52*xp2A2-l/2*B2v*m51A2*m52- 
2/L*B2v*m51 A2*m52*xp2 ; 

L13_3  =  -1/ (2*L) *B2v*mll A2*ml2*xpl-1/ (2*LA3) *B2v*mll A2*ml2*xpl A3; 
L13_4=+l/LA3*B2v*m41*m81*m82+l/ (2*LA3) *B2v*m42*m81A2- 
1/ (2*LA3) *B2v*m41A2*m82 ; 

L13_5=-l/LA3*B2v*m41*m42*m81-2/ (LA2) *B2v*m81*m51*m52*xp2- 
l/L*B2v*m81*m51*m52 ; 

L13_6  =  -1/ (LA3) *B2v*m81*m51*m52*xp2 A2-l/ (LA2) *B2v*m82*m51 A2*xp2; 

L13_7  =  -1/ (2*L) *B2v*m82*m51A2-l/ (2*LA3) *B2v*m82*m51 A2*xp2 A2; 

L13_8=2/ (LA2) *B2v*m41*m51*m52*xp2+l/L*B2v*m41*m51*m52+l/LA3*B2v*m41*m51* 
m52*xp2 A2 ; 

L13_9=l/LA2*B2v*m42*m51 A2*xp2+1/ (2*L) *B2v*m42*m51 A2+l/ (2*LA3) *B2v*m42*m5 
1 A2*xp2 A2 ; 

L13_10  =  -l/LA3*B2v*xpl A2*m81*mll*ml2-1/ (2*LA3) *B2v*xplA2*m82*mllA2; 
L13_ll  =  -l/LA3*B2v*xplA2*m81*mll*ml2-l/ (2*LA3) *B2v*xplA2*m82*mllA2; 
L13_12  =  -1/ (2*LA3) *B2v*xpl*m81A2*ml2-l/LA3*B2v*xpl*m81*m82*mll; 
L13_13=-l/ (2*LA2 ) *B2v*m41 A2*m52-1/ (2*LA2) *B2v*m41A2*m52*xp2- 
l/LA2*B2v*m41*m42*m51 ; 

L13_14  =  -l/LA3*B2v*m41*m42*m51*xp2-l/ (2*LA3) *B2v*xpl*m41 A2*ml2; 

L13_15  =  -l/LA3*B2v*xpl*m41*m42*mll+l/LA3*B2v*xpl A2*m41*mll*ml2 ; 

L13_16  =  +1/ (2*LA3) *B2v*xpl A2*m42*mll A2+l/LA2*B2v*xp2*m41*m81*m52 ; 
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L13_17=l/LA2*B2v*m41*m81*m52+l/LA3*B2v*xp2*m41*m82*m51+l/LA2*B2v*m41*m82 
*m5 1  ; 

L13_18=+l/LA3*B2v*xp2*m42*m81*m51+l/LA2*B2v*m42*m81*m51+l/LA3*B2v*m41*m8 
l*ml2*xpl ; 

L13_19  =+l/LA3*B2v*m41*m82*mll*xpl+l/LA3*B2v*m42*m81*mll*xpl; 

L13_20  =  -l/LA3*B2v*m81*mll*m52*xp2*xpl-l/LA2*B2v*m81*mll*m52*xpl ; 

L13_21  =  -l/LA3*B2v*m81*ml2*m51*xp2*xpl-l/LA2*B2v*m81*ml2*m51*xpl; 

L13_22  =  -l/LA3*B2v*m82*mll*m51*xp2*xpl-l/LA2*B2v*m82*mll*m51*xpl ; 

L13_23  =  +l/LA3*B2v*m41*mll*m52*xp2*xpl+l/LA2*B2v*m41*mll*m52*xpl ; 

L13_24  =  +l/LA3*B2v*m41*ml2*m51*xp2*xpl+l/LA2*B2v*m41*ml2*m51*xpl ; 

L13_25  =  +l/LA3*B2v*m42*mll*m51*xp2*xpl+l/LA2*B2v*m42*mll*m51*xpl ; 

L13_26  =  -l/LA3*B2v*mll*m51*m52*xp2A2*xpl-l/L*B2v*mll*m51*m52*xpl; 
L13_27=-2/LA2*B2v*mll*m51*m52*xp2*xpl-l/ (2*LA3) *B2v*ml2*m51 A2*xp2 A2*xpl; 
L13_28  =  -1/ (2*L) *B2v*ml2*m51A2*xpl-l/LA2*B2v*ml2*m51A2*xp2*xpl; 
L13_29=-l/ (2*LA3) *B2v*mll A2*m52*xp2*xpl A2-l/ (2*LA2) *B2v*mllA2*m52*xplA2; 
L13_30  =  -l/LA3*B2v*mll*ml2*m51*xp2*xpl A2-l/LA2*B2v*mll*ml2*m51*xpl A2 ; 

L13  =  L13_l+L13_2+L13_3+L13_4+L13_5+L13_6+L13_7+L13_8+L13_9+L13_10 . . . 

+L13_11+L13_12+L13_13+L13_14+L13_15+L13_16+L13_17+L13_18+L13_190 
+L13_20+L13_21+L13_22+L13_23+L13_24+L13_25+L13_26+L13_27+L13_28O 
+L13_29+L13_30; 

L14_l=l/2*B2v/LA3*m41*m42A2-l/2*B2v/LA3*m81*m82A2- 
1/ (2*LA3) *B2v*m51*m52 A2*xp2 A3 ; 

L14_2=-3/ (2*LA2 ) *B2v*m51*m52A2*xp2A2-l/2*B2v*m51*m52A2- 
2/L*B2v*m51*m52 A2*xp2 ; 

L14_3=-l/ (2*L) *B2v*mll*ml2 A2*xpl-1/ (2*LA3) *B2v*mll*ml2A2*xplA3; 

L14_4=+l/ (2*LA3) *B2v*m41*m82A2+l/LA3*B2v*m42*m81*m82- 
l/LA3*B2v*m41*m42*m82 ; 

L14_5=-l/ (2*LA3) *B2v*m42 A2*m81-1/ (LA2) *B2v*m81*m52 A2*xp2- 
1/ (2*L) *B2v*m81*m52A2 ; 

L14_6=-l/ (2*LA3) *B2v*m81*m52 A2*xp2 A2-2/ (LA2) *B2v*m82*m51*m52*xp2; 
L14_7=-l/L*B2v*m82*m51*m52- 

l/LA3*B2v*m82*m51*m52*xp2 A2+l/LA2*B2v*m41*m52 A2*xp2 ; 

L14_8=+l/ (2*L) *B2v*m41*m52A2*+l/ (2*LA3) *B2v*m41*m52A2*xp2A2; 

L14_9=+2/ (LA2) *B2v*m42*m51*m52*xp2+l/L*B2v*m42*m51*m52+l/LA3*B2v*m42*m51 
*m52*xp2 A2 ; 

L14_10=— 1/ (2*LA3) *B2v*xplA2*m81*ml2A2-l/LA3*B2v*xplA2*m82*mll*ml2; 

LI 4_1 1=-1 / (2*LA3) *B2v*xpl A2*m81*ml2 A2-l/LA3*B2v*xpl A2*m82*mll*ml2 ; 
L14_12=-l/LA3*B2v*xpl*m81*m82*ml2-l/ (2*LA3) *B2v*xpl*m82 A2*mll; 
L14_13=-l/LA2*B2v*m41*m42*m52-l/LA3*B2v*m41*m42*m52*xp2- 
1/ (2*LA2 ) *B2v*m42 A2*m51 ; 

LI 4_1 4=— 1 / (2*LA3) *B2v*m42 A2*m51*xp2-l/LA3*B2v*xpl*m41*m42*ml2 ; 

LI 4_15=-1 / (2*LA3) *B2v*xpl*m42 A2*mll+1/ (2*LA3) *B2v*xplA2*m41*ml2A2; 
L14_16=+l/LA3*B2v*xpl A2*m42*mll*ml2+l/LA3*B2v*xp2*m41*m82*m52+l/LA2*B2v* 
m41*m82*m52 ; 

L14_17=+l/LA3*B2v*xp2*m42*m81*m52+l/LA2*B2v*m42*m81*m52+l/LA3*B2v*xp2*m4 

2*m82*m51; 

L14_18=+l/LA2*B2v*m42*m82*m51+l/LA3*B2v*m41*m82*ml2*xpl+l/LA3*B2v*m42*m8 
l*ml2*xpl ; 

L14_19=+l/LA3*B2v*m42*m82*mll*xpl-l/LA3*B2v*m81*ml2*m52*xp2*xpl ; 
L14_20=-l/LA2*B2v*m81*ml2*m52*xpl-l/LA3*B2v*m82*mll*m52*xp2*xpl ; 
L14_21=-l/LA2*B2v*m82*mll*m52*xpl-l/LA3*B2v*m82*ml2*m51*xp2*xpl ; 
L14_22=-l/LA2*B2v*m82*ml2*m51*xpl+l/LA3*B2v*m41*ml2*m52*xp2*xpl ; 
L14_23=+l/LA2*B2v*m41*ml2*m52*xpl+l/LA3*B2v*m42*mll*m52*xp2*xpl ; 
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L14_24=+l/LA2*B2v*m42*mll*m52*xpl+l/LA3*B2v*m42*ml2*m51*xp2*xpl ; 
L14_25=+l/LA2*B2v*m42*ml2*m51*xpl-l/ (2*LA3) *B2v*mll*m52A2*xp2A2*xpl; 
L14_26=-l/ (2*L) *B2v*mll*m52 A2*xpl-l/LA2*B2v*mll*m52 A2*xp2*xpl ; 

LI 4_2  7=— 1 / (LA3) *B2v*ml2*m51*m52*xp2 A2*xpl-l/L*B2v*ml2*m51*m52*xpl ; 
L14_28=-2/LA2*B2v*ml2*m51*m52*xp2*xpl-l/ (LA3) *B2v*mll*ml2*m52*xp2*xplA2; 
L14_29=-l/ (LA2) *B2v*mll*ml2*m52*xpl A2-l/ (2*LA3) *B2v*ml2 A2*m51*xp2*xpl A2; 
L14_30=  -l/LA2*B2v*ml2A2*m51*xplA2; 

L14=L14_1+L14_2+L14_3+L14_4+L14_5+L14_6+L14_7+L14_8+L14_9+L14_10 . . . 

+L1 4_1 1+L1 4_12+L1 4_1 3+L1 4_1 4+L1 4_1 5+L1 4_1 6+L1 4_1 7  +L1 4_1 8+L1 4_1 9+L1 4_2  0 .  . 

+L14_21+L14_22+L14_23+L14_24+L14_25+L14_26+L14_27+L14_28+L14_29+L14_30; 

L21_l=l/ (6*LA3) *B2r*m41A3-l/ (6*LA3) *B2r*m81 A3-l/ (6*LA3) *B2r*m51A3*xp2A3; 
L21_2=-l/ (2*LA2 ) *B2r*m51A3*xp2A2-l/6*B2r*m51A3-2/ (3*L) *B2r*m51 A3*xp2; 
L21_3=-l/ (6*L) *B2r*mllA3*xpl- 

1/ (6*LA3) *B2r*mllA3*xplA3+l/ (2*LA3) *B2r*m41*m81A2; 

L2 1_4=-1 / (2*LA3) *B2r*m41A2*m81-l/ (LA2) *B2r*m81*m51A2*xp2- 
1/ (2*L) *B2r*m81*m51A2 ; 

L21_5=-l/ (2*LA3) *B2r*m81*m51 A2*xp2 A2+l/ (LA2) *B2r*m41*m51A2*xp2; 

L21_6=+l/ (2*L) *B2r*m41*m51 A2+l/ (2*LA3) *B2r*m41*m51A2*xp2A2; 

L21_7=-l/ (2*LA3) *B2r*xpl A2*m81*mll A2-l/ (2*LA3) *B2r*xplA2*m81*mllA2; 
L21_8=-l/ (2*LA3) *B2r*xpl*m81 A2*mll-1/ (2*LA2) *B2r*m41 A2*m51; 

L21_9=-l/ (2*LA3) *B2r*m41 A2*m51*xp2-1/ (2*LA3) *B2r*xpl*m41A2*mll; 
L21_10=+l/ (2*LA3) *B2r*xpl A2*m41*mll A2+l/ (LA3) *B2r*xp2*m41*m81*m51; 
L21_ll=+1/ (LA2) *B2r*m41*m81*m51+l/ (LA3) *B2r*m41*m81*mll*xpl; 

L21_12=-l/ (LA3) *B2r*m81*mll*m51*xp2*xpl-l/ (LA2) *B2r*m81*mll*m51*xpl; 
L21_13=+l/ (LA3) *B2r*m41*mll*m51*xp2*xpl+l/ (LA2) *B2r*m41*mll*m51*xpl; 

L2 1_1 4=— 1 / (2*LA3) *B2r*mll*m51A2*xp2A2*xpl-l/ (2*L) *B2r*mll*m51A2*xpl; 
L21_15=-l/ (LA2) *B2r*mll*m51 A2*xp2*xpl-1/ (2*LA3) *B2r*mllA2*m51*xp2*xplA2; 
L21_16=-l/ (2*LA2 ) *B2r*mllA2*m51*xplA2  ; 

L21=  L21_l+L21_2+L21_3+L21_4+L21_5+L21_6+L21_7+L21_8+L21_9+L21_10 . . . 

+L2 1_1 1 +L2 1_1 2  +L2 1_1 3+L2 1_1 4  +L2 1_1 5+L2 1_1 6 ; 

L22_l=l/ (6*LA3) *B2r*m42 A3-l/ (6*LA3) *B2r*m82 A3-l/ (6*LA3) *B2r*m52A3*xp2A3; 
L22_2=-l/ (2*LA2 ) *B2r*m52A3*xp2A2-l/6*B2r*m52A3-2/ (3*L) *B2r*m52 A3*xp2; 
L22_3=-l/ (6*L) *B2r*ml2A3*xpl- 

1/ (6*LA3) *B2r*ml2A3*xplA3+l/ (2*LA3) *B2r*m42*m82 A2; 

L22_4=-l/ (2*LA3) *B2r*m42 A2*m82-1/ (LA2) *B2r*m82*m52 A2*xp2- 
1/ (2*L) *B2r*m82*m52A2; 

L22_5=-l/ (2*LA3) *B2r*m82*m52 A2*xp2 A2+l/ (LA2) *B2r*m42*m52 A2*xp2; 

L22_6=+l/ (2*L) *B2r*m42*m52 A2+l/ (2*LA3) *B2r*m42*m52 A2*xp2 A2; 

L22_7=-l/ (2*LA3) *B2r*xpl A2*m82*ml2 A2-l/ (2*LA3) *B2r*xpl A2*m82*ml2 A2; 
L22_8=-l/ (2*LA3) *B2r*xpl*m82 A2*ml2-1/ (2*LA2) *B2r*m42 A2*m52; 

L22_9=-l/ (2*LA3) *B2r*m42 A2*m52*xp2-1/ (2*LA3) *B2r*xpl*m42 A2*ml2; 
L22_10=+l/ (2*LA3) *B2r*xpl A2*m42*ml2 A2+l/ (LA3) *B2r*xp2*m42*m82*m52; 
L22_ll=+1/ (LA2) *B2r*m42*m82*m52+l/ (LA3) *B2r*m42*m82*ml2*xpl; 

L22_12=-l/ (LA3) *B2r*m82*ml2*m52*xp2*xpl-l/ (LA2) *B2r*m82*ml2*m52*xpl; 
L22_13=  +1/ (LA3) *B2r*m42*ml2*m52*xp2*xpl+l/ (LA2) *B2r*m42*ml2*m52*xpl; 
L22_14=-l/ (2*LA3) *B2r*ml2*m52 A2*xp2 A2*xpl-1/ (2*L) *B2r*ml2*m52 A2*xpl; 
L22_15  =-1/ (LA2) *B2r*ml2*m52 A2*xp2*xpl- 
1/ (2*LA3) *B2r*ml2 A2*m52*xp2*xpl A2 ; 

L22_l 6  =  -1/ (2*LA2) *B2r*ml2 A2*m52*xpl A2; 
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L22 


=  L22_l+L22_2+L22_3+L22_4+L22_5+L22_6+L22_7+L22_8+L22_9+L22_10 . . . 
+L22_11+L22_12+L22_13+L22_14+L22_15+L22_16; 

L23_l=l/ (2*LA3) *B2r*m41A2*m42-l/ (2*LA3) *B2r*m81 A2*m82- 
1/ (2*LA3) *B2r*m51A2*m52*xp2A3; 

L23_2=-3/ (2*LA2 ) *B2r*m51 A2*m52*xp2 A2-l/2*B2r*m51 A2*m52- 
2/L*B2r*m51 A2*m52*xp2; 

L23_3=-l/ (2*L) *B2r*mllA2*ml2*xpl-l/ (2*LA3) *B2r*mllA2*ml2*xplA3; 

L23_4=+l/ (LA3) *B2r*m41*m81*m82+l/ (2*LA3) *B2r*m42*m81 A2- 
1/ (2*LA3) *B2r*m41A2*m82  ; 

L23_5=-l/ (LA3) *B2r*m41*m42*m81-2/ (LA2) *B2r*m81*m51*m52*xp2- 
l/L*B2r*m81*m51*m52 ; 

L23_6=-l/ (LA3) *B2r*m81*m51*m52*xp2 A2-l/ (2*LA2) *B2r*m82*m51 A2*xp2; 
L23_7=-l/ (2*L) *B2r*m82*m51 A2-l/ (2*LA3) *B2r*m82*m51 A2*xp2 A2; 

L23_8=+2/ (LA2) *B2r*m41*m51*m52*xp2+l/L*B2r*m41*m51*m52; 

L23_9=+l/ (LA3) *B2r*m41*m51*m52*xp2 A2+l/ (LA2) *B2r*m42*m51 A2*xp2; 
L23_10=+l/ (2*L) *B2r*m42*m51 A2+l/ (2*LA3) *B2r*m42*m51 A2*xp2 A2; 

L23_ll=-1/ (LA3) *B2r*xplA2*m81*mll*ml2-l/ (2*LA3) *B2r*xpl A2*m82*mll A2; 
L23_12=-l/ (LA3) *B2r*xpl A2*m81*mll*ml2-1/ (2*LA3) *B2r*xpl A2*m82*mll A2; 
L23_13=-l/ (2*LA3) *B2r*xpl*m81 A2*ml2-1/ (LA3) *B2r*xpl*m81*m82*mll; 
L23_14=-l/ (2*LA2 ) *B2r*m41A2*m52-l/ (2*LA3) *B2r*m41A2*m52*xp2; 

L23_15=-l/ (LA2) *B2r*m41*m42*m51-l/ (LA3) *B2r*m41*m42*m51*xp2; 

L23_16=-l/ (2*LA3) *B2r*xpl*m41A2*ml2-l/ (LA3) *B2r*xpl*m41*m42*mll; 
L23_17=+l/ (LA3) *B2r*xpl A2*m41*mll*ml2+1/ (2*LA3) *B2r*xpl A2*m42*mll A2; 
L23_18=+l/ (LA3) *B2r*xp2*m41*m81*m52+l/ (LA2) *B2r*m41*m81*m52; 

L23_19=+l/ (LA3) *B2r*xp2*m41*m82*m51+l/ (LA2) *B2r*m41*m82*m51; 

L23_20=+l/ (LA3) *B2r*xp2*m42*m81*m51+l/ (LA2) *B2r*m42*m81*m51; 

L23_21=+l/ (LA3) *B2r*m41*m81*ml2*xpl+l/ (LA3) *B2r*m41*m82*mll*xpl; 
L23_22=+l/ (LA3) *B2r*m42*m81*mll*xpl-l/ (LA3) *B2r*m81*mll*m52*xp2*xpl; 
L23_23=-l/ (LA2) *B2r*m81*mll*m52*xpl-l/ (LA3) *B2r*m81*ml2*m51*xp2*xpl; 
L23_24=-l/ (LA2) *B2r*m81*ml2*m51*xpl-l/ (LA3) *B2r*m82*mll*m51*xp2*xpl; 
L23_25=-l/ (LA2) *B2r*m82*mll*m51*xpl+l/ (LA3) *B2r*m41*mll*m52*xp2*xpl; 
L23_26=+l/ (LA2) *B2r*m41*mll*m52*xpl+l/ (LA3) *B2r*m41*ml2*m51*xp2*xpl; 
L23_27=+l/ (LA2) *B2r*m41*ml2*m51*xpl+l/ (LA3) *B2r*m42*mll*m51*xp2*xpl; 
L23_28=+l/ (LA2) *B2r*m42*mll*m51*xpl-l/ (LA3) *B2r*mll*m51*m52*xp2A2*xpl; 
L23_29=-l/L*B2r*mll*m51*m52*xpl-2/ (LA2) *B2r*mll*m51*m52*xp2*xpl; 
L23_30=-l/ (2*LA3) *B2r*ml2*m51 A2*xp2 A2*xpl-1/ (2*L) *B2r*ml2*m51 A2*xpl; 
L23_31=-l/ (LA2) *B2r*ml2*m51 A2*xp2*xpl-1/ (2*LA3) *B2r*mllA2*m52*xp2*xplA2; 
L23_32=-l/ (2*LA2 ) *B2r*mllA2*m52*xplA2-l/ (LA3) *B2r*mll*ml2*m51*xp2*xplA2; 
L23_33=-l/ (LA2) *B2r*mll*ml2*m51*xplA2; 


L23=L23_1+L23_2+L23_3+L23_4+L23_5+L23_6+L23_7+L23_8+L23_9+L23_10 . . . 

+L23_11+L23_12+L23_13+L23_14+L23_15+L23_16+L23_17+L23_18+L23_19+L23_20 . . 

+L23_21+L23_22+L23_23+L23_24+L23_25+L23_26+L23_27+L23_28+L23_29+L23_30 . . 
+L2 3_3 1 +L2  3_3  2  +L2  3_3  3 ; 


L2  4_1  =  1 / (2*LA3) *B2r*m41*m42 A2-l/ (2*LA3) *B2r*m81*m82 A2- 
1/ (2*LA3) *B2r*m51*m52 A2*xp2 A3; 


43 


L24_2=-3/ (2*LA2 ) *B2r*m51*m52A2*xp2A2-l/2*B2r*m51*m52A2- 
2/L*B2r*m51*m52 A2*xp2 ; 

L24_3  =  -1/ (2*L) *B2r*mll*ml2A2*xpl-l/ (2*LA3) *B2r*mll*ml2A2*xplA3; 
L24_4=+l/ (2*LA3) *B2r*m41*m82A2+l/ (LA3) *B2r*m42*m81*m82- 
1/ (LA3) *B2r*m41*m42*m82 ; 

L24_5=-l/ (2*LA3) *B2r*m42A2*m81-l/ (LA2) *B2r*m81*m52A2*xp2- 
1/ (2*L) *B2r*m81*m52A2 ; 

L24_6=-l/ (2*LA3) *B2r*m81*m52A2*xp2A2-l/ (2*LA2) *B2r*m82*m51*m52*xp2; 
L24_7=-l/L*B2r*m82*m51*m52- 

1/ (LA3) *B2r*m82*m51*m52*xp2 A2+l/ (LA2) *B2r*m41*m52A2*xp2; 

L24_8=+l/ (2*L) *B2r*m41*m52A2+l/ (2*LA3) *B2r*m41*m52A2*xp2A2; 

L24_9=+2/ (LA2) *B2r*m42*m51*m52*xp2+l/L*B2r*m42*m51*m52+l/ (LA3) *B2r*m42*m 
51*m52*xp2A2; 

L24_10=-l/ (2*LA3) *B2r*xpl A2*m81*ml2 A2-l/ (LA3) *B2r*xplA2*m82*mll*ml2; 

L2  4_1 1=-1 / (2*LA3) *B2r*xpl A2*m81*ml2 A2-l/ (LA3) *B2r*xplA2*m82*mll*ml2; 
L24_12=-l/ (LA3) *B2r*xpl*m81*m82*ml2-l/ (2*LA3) *B2r*xpl*m82 A2*mll; 
L24_13=-l/ (LA2) *B2r*m41*m42*m52-l/ (LA3) *B2r*m41*m42*m52*xp2- 
1/ (2*LA2 ) *B2r*m42 A2*m51 ; 

L2  4_1 4=— 1 / (2*LA3) *B2r*m42 A2*m51*xp2-1/ (LA3) *B2r*xpl*m41*m42*ml2; 
L24_15=-l/ (2*LA3) *B2r*xpl*m42A2*mll+l/ (2*LA3) *B2r*xplA2*m41*ml2A2; 
L24_16=+l/ (LA3) *B2r*xpl A2*m42*mll*ml2+1/ (LA3) *B2r*xp2*m41*m82*m52; 
L24_17=+l/ (LA2) *B2r*m41*m82*m52+l/ (LA3) *B2r*xp2*m42*m81*m52+l/ (LA2) *B2r* 
m42*m81*m52 ; 

L24_18=+l/ (LA3) *B2r*xp2*m42*m82*m51+l/ (LA2) *B2r*m42*m82*m51; 

L24_19=+l/ (LA3) *B2r*m41*m82*ml2*xpl+l/ (LA3) *B2r*m42*m81*ml2*xpl; 
L24_20=+l/ (LA3) *B2r*m42*m82*mll*xpl-l/ (LA3) *B2r*m81*ml2*m52*xp2*xpl; 
L24_21=-l/ (LA2) *B2r*m81*ml2*m52*xpl-l/ (LA3) *B2r*m82*mll*m52*xp2*xpl; 
L24_22=-l/ (LA2) *B2r*m82*mll*m52*xpl-l/ (LA3) *B2r*m82*ml2*m51*xp2*xpl; 
L24_23=-l/ (LA2) *B2r*m82*ml2*m51*xpl+l/ (LA3) *B2r*m41*ml2*m52*xp2*xpl; 
L24_24=+l/ (LA2) *B2r *m4 1 *ml2 *m52 *xpl+l / (LA3) *B2r*m42*mll*m52*xp2*xpl; 
L24_25=+l/ (LA2) *B2r*m42*mll*m52*xpl+l/ (LA3) *B2r*m42*ml2*m51*xp2*xpl; 
L24_26=+l/ (LA2) *B2r*m42*ml2*m51*xpl-l/ (2*LA3) *B2r*mll*m52A2*xp2A2*xpl; 

L2  4_2  7=— 1 / (2*L) *B2r*mll*m52A2*xpl-l/ (LA2) *B2r*mll*m52A2*xp2*xpl; 
L24_28=-l/ (LA3) *B2r*ml2*m51*m52*xp2 A2*xpl-l/L*B2r*ml2*m51*m52*xpl; 
L24_29=-2/ (LA2) *B2r*ml2*m51*m52*xp2*xpl- 
1/ (LA3) *B2r*mll*ml2*m52*xp2*xpl A2 ; 

L24_30=-l/ (LA2) *B2r*mll*ml2*m52*xpl A2-l/ (2*LA3) *B2r*ml2 A2*m51*xp2*xpl A2; 
L24_31=-l/ (2*LA2 ) *B2r*ml2 A2*m51*xpl A2 ; 

L2  4=  L24_l+L24_2+L24_3+L24_4+L24_5+L24_6+L24_7+L24_8+L24_9+L24_10 .  .  . 

+L24_11+L24_12+L24_13+L24_14+L24_15+L24_16+L24_17+L24_18+L24_19+L24_20 . . 

+L24_21+L24_22+L24_23+L24_24+L24_25+L24_26+L24_27+L24_28+L24_29+L24_30 . . 
+L24_31; 

L31_l=l/ (6*LA3) *Blv*m41A3-l/ (6*LA3) *Blv*m81 A3-l/ (6*L) *Blv*m51 A3*xp2; 
L31_2=-l/ (6*LA3) *Blv*m51A3*xp2A3-l/ (2*LA2) *Blv*mllA3*xplA2- 
1/ 6*Blv*mllA3; 

L31_3=-2/ (3*L) *Blv*mllA3*xpl- 

1/ (6*LA3) *Blv*mll A3*xpl A3+l/ (2*LA3) *Blv*m41*m81A2; 

L31_4=+l/ (2*LA3) *Blv*m41*m81A2-l/ (2*LA3) *Blv*xp2A2*m81*m51A2; 

L31_5=+l/ (2*LA3) *Blv*xp2 A2*m41*m51 A2-l/ (LA2) *Blv*m81*mll A2*xpl; 

L31_6=-l/ (2*LA3) *Blv*m81*mll A2*xpl A2-l/ (2*L) *Blv*m81*mllA2; 

L31_7=-l/ (2*LA3) *Blv*xp2*m81 A2*m51-1/ (2*LA2) *Blv*m81A2*mll; 
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L31_8=-l/ (2*LA3) *Blv*m81A2*mll*xpl-l/ (2*LA3) *Blv*xp2*m41A2*m51; 

L31_9=-l/ (2*LA2 ) *Blv*m41A2*mll- 

1/ (2*LA3) *Blv*m41A2*mll*xpl+l/ (2*L) *Blv*m41*mllA2; 

L31_10=+l/ (LA2) *Blv*m41*mllA2*xpl+l/ (2*LA3) *Blv*m41*mllA2*xplA2; 
L31_ll=+1/ (LA3) *Blv*xp2*m41*m81*m51+l/ (LA2) *Blv*m41*m81*mll; 

L31_12=+l/ (LA3) *Blv*m41*m81*mll*xpl-l/ (LA3) *Blv*m81*mll*m51*xp2*xpl; 
L31_13=-l/ (LA2) *Blv*m81*mll*m51*xp2+l/ (LA3) *Blv*m41*mll*m51*xp2*xpl; 
L31_14=+l/ (LA2) *Blv*m41*mll*m51*xp2-l/ (2*LA2) *Blv*mll*m51A2*xp2/'2; 
L31_15=-l/ (2*LA3) *Blv*mll*m51 A2*xp2A2*xpl-l/ (2*L) *Blv*mll A2*m51*xp2; 
L31_16=-l/ (LA2) *Blv*mllA2*m51*xp2*xpl-l/ (2*LA3) *Blv*mllA2*m51*xp2*xpl"2; 
L31_17— -1/ 3*Clv/ deltas at A2*kpsiA3*mllA3-l/3*Clv/ delta s at A2*kv~3*m2 1A3; 
L31_18=-l/3*Clv/deltasatA2*kr/''3*m31A3-l/3*Clv/deltasatA2*kyA3*m41A3; 
L31_l 9=-Clv/ delta sat A2*kpsi A2*kv*m2 l*ml 1 A2- 
Clv/ delta sat A2*kpsi A2*kr*m31*ml 1 A2 ; 

L31_20=-Clv/ delta sat A2*kpsiA2*ky*m41*mll A2- 
Clv/ delta sat A2*kpsi*kvA2*m21 A2*ml 1 ; 

L31_21=-Clv/ deltasat A2*kpsi*kr A2*m31 A2*mll- 
Civ/ delta sat A2*kpsi*kyA2*m4 1 A2*ml 1 ; 
L31_22=-Clv/deltasatA2*kvA2*kr*m31*m21A2- 
Clv/ delta sat A2*kvA2*ky*m4 l*m21 A2 ; 
L31_23=-Clv/deltasatA2*kv*krA2*m21*m31A2- 
Clv/ delta sat A2*kv*kyA2*m21*m4 1 A2 ; 

L31_24  =  -Civ/ deltasat A2*kr A2*ky*m41*m31A2- 
Clv/ delta sat A2*kr*kyA2*m31*m4 1 A2 ; 

L31_25  =  -2*Clv/deltasat A2*kpsi*kv*kr*m31 *m21 *ml 1 ; 

L31_26  =  -2*Clv/deltasat A2*kpsi*kv*ky*m4 l*m2 l*ml 1 ; 

L31_27  =  -2*Clv/deltasat A2*kpsi*kr*ky*m4 l*m31 *ml 1 ; 

L31_28  =  -2*Clv/ deltasat A2*kv*kr*ky*m4 l*m2 l*m31 ; 

L31  =  L31_l+L31_2+L31_3+L31_4+L31_5+L31_6+L31_7+L31_8+L31_9+L31_10 . . . 

+L31_11+L31_12+L31_13+L31_14+L31_15+L31_16+L31_17+L31_18+L31_19+L31_20 . . 
+L31_21+L31_22+L31_23+L31_24+L31_25+L31_26+L31_27+L31_28; 


L32_l=l/ (6*LA3) *Blv*m42A3-l/ (6*LA3) *Blv*m82A3-l/ (6*L) *Blv*m52 A3*xp2; 
L32_2=-l/ (6*LA3) *Blv*m52 A3*xp2 A3-l/ (2*LA2) *Blv*ml2 A3*xpl A2- 
1/ 6*Blv*ml2A3; 

L32_3=-2/ (3*L) *Blv*ml2 A3*xpl- 

1/ (6*LA3) *Blv*ml2 A3*xpl A3+l/ (2*LA3) *Blv*m42*m82A2; 

L32_4=+l/ (2*LA3) *Blv*m42*m82 A2-l/ (2*LA3) *Blv*xp2 A2*m82*m52 A2; 

L32_5=+l/ (2*LA3) *Blv*xp2 A2*m42*m52 A2-l/ (LA2) *Blv*m82*ml2 A2*xpl; 

L32_6=-l/ (2*LA3) *Blv*m82*ml2 A2*xpl A2-l/ (2*L) *Blv*m82*ml2 A2; 

L32_7=-l/ (2*LA3) *Blv*xp2*m82 A2*m52-1/ (2*LA2) *Blv*m82 A2*ml2; 

L32_8=-l/ (2*LA3) *Blv*m82 A2*ml2*xpl-1/ (2*LA3) *Blv*xp2*m42 A2*m52; 

L32_9=-l/ (2*LA2 ) *Blv*m42 A2*ml2- 

1/ (2*LA3) *Blv*m42A2*ml2*xpl+l/ (2*L) *Blv*m42*ml2A2; 

L32_10=+l/ (LA2) *Blv*m42*ml2 A2*xpl+1/ (2*LA3) *Blv*m42*ml2 A2*xpl A2; 
L32_ll=+1/ (LA3) *Blv*xp2*m42*m82*m52+l/ (LA2) *Blv*m42*m82*ml2; 

L32_12=+l/ (LA3) *Blv*m42*m82*ml2*xpl-l/ (LA3) *Blv*m82*ml2*m52*xp2*xpl; 
L32_13=-l/ (LA2) *Blv*m82*ml2*m52*xp2+l/ (LA3) *Blv*m42*ml2*m52*xp2*xpl; 
L32_14=+l/ (LA2) *Blv*m42*ml2*m52*xp2-l/ (2*LA2) *Blv*ml2*m52A2*xp2 A2; 
L32_15=-l/ (2*LA3) *Blv*ml2*m52 A2*xp2 A2*xpl-1/ (2*L) *Blv*ml2 A2*m52*xp2; 
L32_16=-l/ (LA2) *Blv*ml2 A2*m52*xp2*xpl-1/ (2*LA3) *Blv*ml2 A2*m52*xp2*xpl A2; 
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L32_17=-l/ 3*Clv/ deltas at A2*kpsiA3*ml2A3-l/3*Clv/ deltasat A2*kvA3*m22 A3; 
L32_18=-l/3*Clv/ deltasat A2*kr A3*m32A3-l/3*Clv/ deltasat A2*kyA3*m42 A3; 
L32_l 9=-Clv/ delta sat A2*kpsi A2*kv*m22*ml2A2- 
Clv/ delta sat A2*kpsiA2*kr*m32*ml2 A2 ; 

L32_20=-Clv/ delta sat A2*kpsiA2*ky*m42*ml2 A2- 
Cl v/ deltasatA2*kpsi*kvA2*m22 A2*ml2 ; 
L32_21=-Clv/deltasatA2*kpsi*krA2*m32A2*ml2- 
Civ/ delta sat A2*kpsi*kyA2*m42 A2*ml2 ; 

L32_22=-Clv/ deltasat A2*kvA2*kr*m32*m22 A2- 
Clv/ delta sat A2*kvA2*ky*m42*m22 A2 ; 

L32_23=-Clv/ deltasat A2*kv*kr A2*m22*m32 A2- 
Clv/ delta sat A2*kv*kyA2*m22*m42 A2 ; 

L32_24=-Clv/ deltasat A2*kr A2*ky*m42*m32 A2- 
Clv/ delta sat A2*kr*kyA2*m32*m42 A2 ; 

L32_25  =  -2*Clv/deltasat A2*kpsi*kv*kr*m32*m22*ml2; 

L32_26  =  -2*Clv/deltasat A2*kpsi*kv*ky*m42*m22*ml2; 

L32_27  =  -2*Clv/deltasat A2*kpsi*kr*ky*m42*m32*ml2; 

L32_28  =  -2*Clv/deltasat A2*kv*kr*ky*m42*m22*m32 ; 


L32  =  L32_l+L32_2+L32_3+L32_4+L32_5+L32_6+L32_7+L32_8+L32_9+L32_10 . . . 
+L32_11+L32_12+L32_13+L32_14+L32_15+L32_16+L32_17+L32_18+L32_19+L32_20 . . 
+L32_21+L32_22+L32_23+L32_24+L32_25+L32_26+L32_27+L32_28; 


L33_l=l/ (2*LA3) *Blv*m41A2*m42-l/ (2*LA3) *Blv*m81A2*m82- 
1/ (2*L) *Blv*m51 A2*m52*xp2 ; 

L33_2=-l/ (2*LA3) *Blv*m51 A2*m52*xp2A3-3/ (2*LA2) *Blv*mllA2*ml2*xplA2; 

L33_3=-l/2*Blv*mllA2*ml2-2/L*Blv*mllA2*ml2*xpl- 

1/ (2*LA3) *Blv*mllA2*ml2*xplA3; 

L33_4=+l/ (LA3) *Blv*m41*m81*m82+l/ (2*LA3) *Blv*m42*m81 A2+l/ (LA3) *Blv*m41*m 
81*m82  ; 

L33_5=+l/ (2*LA3) *Blv*m42*m81A2-l/ (LA3) *Blv*xp2A2*m81*m51*m52; 

L33_6=-l/ (2*LA3) *Blv*xp2 A2*m82*m51 A2  +  l/ (LA3) *Blv*xp2A2*m41*m51*m52; 
L33_7=+l/ (2*LA3) *Blv*xp2 A2*m42*m51 A2-2/ (LA2) *Blv*m81*mll*ml2*xpl; 
L33_8=-l/ (LA3) *Blv*m81*mll*ml2*xplA2-l/L*Blv*m81*mll*ml2- 
1/ (LA2) *Blv*m82*mllA2*xpl  ; 

L33_9=-l/ (2*LA3) *Blv*m82*mll A2*xpl A2-l/ (2*L) *Blv*m82*mll A2; 

L33_10=-l/ (2*LA3) *Blv*xp2*m81 A2*m52-1/ (LA3) *Blv*xp2*m81*m82*m51; 
L33_ll=-1/ (2*LA2 ) *Blv*m81 A2*ml2-1/ (2*LA3) *Blv*m81A2*ml2*xpl; 

L33_12=-l/ (LA2) *Blv*m81*m82*mll-l/ (LA3) *Blv*m81*m82*mll*xpl  ; 

L33_13=-l/ (2*LA3) *Blv*xp2*m41 A2*m52-1/ (LA3) *Blv*xp2*m41*m42*m51; 
L33_14=-l/ (2*LA2 ) *Blv*m41A2*ml2-l/ (2*LA3) *Blv*m41A2*ml2*xpl- 
1/ (LA2) *Blv*m41*m42*mll; 

L33_15=l/ (LA3) *Blv*m41*m42*mll*xpl+l/L*Blv*m41*mll*ml2+2/ (LA2) *Blv*m41*m 
ll*ml2*xpl; 

L33_16=+l/ (LA3) *Blv*m41*mll*ml2*xpl A2+l/ (2*L) *Blv*m42*mllA2; 

L33_17=+l/ (LA2) *Blv*m42*mll A2*xpl+1/ (2*LA3) *Blv*m42*mllA2*xplA2; 
L33_18=+l/ (LA3) *Blv*xp2*m41*m81*m52+l/ (LA3) *Blv*xp2*m41*m82*m51; 
L33_19=+l/ (LA3) *Blv*xp2*m42*m81*m51+l/ (LA2) *Blv*m41*m81*ml2; 

L33_20=+l/ (LA3) *Blv*m41*m81*ml2*xpl+l/ (LA2) *Blv*m41*m82*mll; 

L33_21=+l/ (LA3) *Blv*m41*m82*mll*xpl+l/ (LA2) *Blv*m42*m81*mll; 

L33_22=+l/ (LA3) *Blv*m42*m81*mll*xpl-l/ (LA3) *Blv*m81*mll*m52*xp2*xpl; 
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L33_23=-l/ (LA2) *Blv*m81*mll*m52*xp2-l/ (LA3) *Blv*m81*ml2*m51*xp2*xpl; 
L33_24  =-1/ (LA2) *Blv*m81*ml2*m51*xp2-l/ (LA3) *Blv*m82*mll*m51*xp2*xpl; 
L33_25=-l/ (LA2) *Blv*m82*mll*m51*xp2+l/ (LA3) *Blv*m41*mll*m52*xp2*xpl; 
L33_26=+l/ (LA2) *Blv*m41*mll*m52*xp2+l/ (LA3) *Blv*m41*ml2*m51*xp2*xpl; 
L33_27=+l/ (LA2) *Blv*m41*ml2*m51*xp2+l/ (LA3) *Blv*m42*mll*m51*xp2*xpl; 
L33_28=+l/ (LA2) *Blv*m42*mll*m51*xp2-l/ (LA2) *Blv*mll*m51*m52*xp2A2  ; 
L33_29=-l/ (LA3) *Blv*mll*m51*m52*xp2A2*xpl-l/ (2*LA2) *Blv*ml2*m51 A2*xp2 A2 
L33_30=-l/ (2*LA3) *Blv*ml2*m51 A2*xp2 A2*xpl-1/ (2*L) *Blv*mll A2*m52*xp2; 
L33_31=-l/ (LA2) *Blv*mll A2*m52*xp2*xpl-1/ (2*LA3) *Blv*mllA2*m52*xp2*xplA2 
L33_32=-l/L*Blv*mll*ml2*m51*xp2-2/ (LA2) *Blv*mll*ml2*m51*xp2*xpl; 
L33_33=-l/ (LA3) *Blv*mll*ml2*m51*xp2*xplA2- 
Clv/deltasatA2*kpsiA3*mllA2*ml2; 

L33_34=-Clv/ delta sat A2*kvA3*m21 A2*m22-Clv/ delta sat A2*kr A3*m31 A2*m32  ; 
L33_35=-Clv/ deltasat A2*kyA3*m41 A2*m42- 
2*C1 v/ delta sat A2*kpsi A2*kv*m2 1 *ml l*ml2 ; 

L33_36=-Clv/ delta sat A2*kpsiA2*kv*m22*mll A2- 
2*Clv/ delta sat A2*kpsi A2*kr*m31*ml l*ml2 ; 

L33_37=-Clv/ delta sat A2*kpsi A2*kr*m32*ml 1 A2- 
2*Clv/ delta sat A2*kpsi A2*ky*m4 l*ml l*ml2 ; 

L33_38=-Clv/ delta sat A2*kpsi A2*ky*m42*ml 1 A2- 
Clv/ delta sat A2*kpsi*kvA2*m21 A2*ml2  ; 
L33_39=-2*Clv/deltasatA2*kpsi*kvA2*m21*m22*mll- 
Clv/ delta sat A2*kpsi*kr A2*m31 A2*ml2 ; 

L33_40=-2*Clv/ delta sat A2*kpsi*kr A2*m31*m32*mll- 
Civ/ delta sat A2*kpsi*kyA2*m4 1 A2*ml2 ; 
L33_41=-2*Clv/deltasatA2*kpsi*kyA2*m41*m42*mll- 
2*Clv/ delta sat A2*kvA2*kr*m31*m21*m22; 

L33_42=-Clv/ deltasat A2*kvA2*kr*m32*m21 A2- 
2*Clv/ delta sat A2*kvA2*ky*m4 I*m21*m22 ; 
L33_43=-Clv/deltasatA2*kvA2*ky*m42*m21A2- 
2*Clv/ delta sat A2*kv*kr A2*m2 I*m31*m32 ; 

L33_44=-Clv/ deltasat A2*kv*krA2*m22 *m31A2- 
2*Clv/ delta sat A2*kv*kyA2*m21*m4 l*m42 ; 

L33_45=-Clv/ delta sat A2*kv*kyA2*m22*m41 A2- 
2*C1 v/ delta sat A2*kr A2*ky*m4 I*m31*m32 ; 

L33_4  6=-Clv/ delta sat A2*kr A2*ky*m42*m31 A2- 
2*Clv/ delta sat A2*kr*kyA2*m31*m4 l*m42 ; 

L33_47=-Clv/ deltasat A2*kr*kyA2*m32*m41 A2- 
2*Clv/ delta sat A2*kpsi*kv*kr*m31*m21*ml2 ; 

L33_48=-2*Clv/ delta sat A2*kpsi*kv*kr*m31*m22*ml 1 ; 

L33_4  9=-2*Clv/ delta sat A2*kpsi*kv*kr*m32*m2 l*ml 1 ; 

L33_50=-2*Clv/ delta sat A2*kpsi*kv*ky*m4 l*m2 l*ml2 ; 

L33_51=-2*Clv/ delta sat A2*kpsi*kv*ky*m4 l*m22*ml 1 ; 

L33_52=-2*Clv/ delta sat A2*kpsi*kv*ky*m42*m2 l*ml 1 ; 

L33_53=-2*Clv/ delta sat A2*kpsi*kr*ky*m4 I*m31*ml2 ; 

L33_54=-2*Clv/ delta sat A2*kpsi*kr*ky*m4 l*m32*ml 1 ; 

L33_55=-2*Clv/ delta sat A2*kpsi*kr*ky*m42*m31*ml 1 ; 

L33_56=-2*Clv/ delta sat A2*kv*kr*ky*m41*m21*m32- 
2*Clv/ delta sat A2*kv*kr*ky*m41*m22*m31 ; 
L33_57=-2*Clv/deltasatA2*kv*kr*ky*m42*m21*m31; 


L33  =  L33_l+L33_2+L33_3+L33_4+L33_5+L33_6+L33_7+L33_8+L33_9+L33_10 . . . 

+L33_11+L33_12+L33_13+L33_14+L33_15+L33_16+L33_17+L33_18+L33_19+L33_20 . 
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+L33_21+L33_22+L33_23+L33_24+L33_25+L33_26+L33_27+L33_28+L33_29+L33_30 . . 


+L33_31+L33_32+L33_33+L33_34+L33_35+L33_36+L33_37+L33_38+L33_39+L33_40 . . 

+L33_41+L33_42+L33_44+L33_45+L33_46+L33_47+L33_48+L33_49+L33_50+L33_51 . . 
+L33_52+L33_53+L33_54+L33_55+L33_56+L33_57 ; 


L34_l=l/ (2*LA3) *Blv*m41*m42 A2-l/ (2*LA3) *Blv*m81*m82A2- 
1/ (2*L) *Blv*m51*m52 A2*xp2 ; 

L34_2=-l/ (2*LA3) *Blv*m51*m52 A2*xp2 A3-3/ (2*LA2) *Blv*mll*ml2A2*xplA2; 

L34_3=-l/2*Blv*mll*ml2A2-2/L*Blv*mll*ml2A2*xpl- 

1/ (2*LA3) *Blv*mll*ml2A2*xplA3; 

L34_4=+l/ (2*LA3) *Blv*m41*m82 A2+l/ (LA3) *Blv*m42*m81*m82+l/ (2*LA3) *Blv*m41 
*m82A2; 

L34_5=+l/ (LA3) *Blv*m42*m81*m82-l/ (2*LA3) *Blv*xp2A2*m81*m52A2; 

L34_6=-l/ (LA3) *Blv*xp2 A2*m82*m51*m52+1/ (2*LA3) *Blv*xp2A2*m41*m52A2; 
L34_7=+l/ (LA3) *Blv*xp2 A2*m42*m51*m52-1/ (LA2) *Blv*m81*ml2 A2*xpl; 

L34_8=-l/ (2*LA3) *Blv*m81*ml2A2*xplA2-l/ (2*L) *Blv*m81*ml2A2; 

L34_9=-2/ (LA2) *Blv*m82*mll*ml2*xpl-l/ (LA3) *Blv*m82*mll*ml2*xplA2; 
L34_10=-l/L*Blv*m82*mll*ml2-l/ (LA3) *Blv*xp2*m81*m82*m52- 
1/ (2*LA3) *Blv*xp2*m82 A2*m51 ; 

L34_ll=-1/ (LA2) *Blv*m81*m82*ml2-l/ (LA3) *Blv*m81*m82*ml2*xpl- 
1/ (2*LA2 ) *Blv*m82A2*mll ; 

L34_12=-l/ (2*LA3) *Blv*m82A2*mll*xpl-l/ (LA3) *Blv*xp2*m41*m42*m52; 
L34_13=-l/ (2*LA3) *Blv*xp2*m42 A2*m51-1/ (LA2) *Blv*m41*m42*ml2; 

L34_14=-l/ (LA3) *Blv*m41*m42*ml2*xpl-l/ (2*LA2) *Blv*m42A2*mll; 

L34_15=-l/ (2*LA3) *Blv*m42A2*mll*xpl+l/ (2*L) *Blv*m41*ml2A2+ 

1/ (LA2) *Blv*m41*ml2 A2*xpl ; 

L34_16=+l/ (2*LA3) *Blv*m41*ml2A2*xplA2+l/L*Blv*m42*mll*ml2; 

L34_17=+2/ (LA2) *Blv*m42*mll*ml2*xpl+l/ (LA3) *Blv*m42*mll*ml2*xplA2; 
L34_18=+l/ (LA3) *Blv*xp2*m41*m82*m52+l/ (LA3) *Blv*xp2*m42*m81*m52; 
L34_19=+l/ (LA3) *Blv*xp2*m42*m82*m51+l/ (LA2) *Blv*m41*m82*ml2; 

L34_20=+l/ (LA3) *Blv*m41*m82*ml2*xpl+l/ (LA2) *Blv*m42*m81*ml2; 

L34_21=+l/ (LA3) *Blv*m42*m81*ml2*xpl+l/ (LA2) *Blv*m42*m82*mll; 

L34_22=+l/ (LA3) *Blv*m42*m82*mll*xpl-l/ (LA3) *Blv*m81*ml2*m52*xp2*xpl; 
L34_23=-l/ (LA2) *Blv*m81*ml2*m52*xp2-l/ (LA3) *Blv*m82*mll*m52*xp2*xpl; 
L34_24=-l/ (LA2) *Blv*m82*mll*m52*xp2-l/ (LA3) *Blv*m82*ml2*m51*xp2*xpl; 
L34_25=-l/ (LA2) *Blv*m82*ml2*m51*xp2+l/ (LA3) *Blv*m41*ml2*m52*xp2*xpl; 
L34_2  6=  +1/ (LA2) *Blv*m41*ml2*m52*xp2+l/ (LA3) *Blv*m42*mll*m52*xp2*xpl; 
L34_27=+l/ (LA2) *Blv*m42*mll*m52*xp2+l/ (LA3) *Blv*m42*ml2*m51*xp2*xpl; 
L34_28=+l/ (LA2) *Blv*m42*ml2*m51*xp2-l/ (2*LA2) *Blv*mll*m52A2*xp2A2; 
L34_29=-l/ (2*LA3) *Blv*mll*m52A2*xp2A2*xpl-l/ (LA2) *Blv*ml2*m51*m52*xp2 A2; 
L34_30=-l/ (LA3) *Blv*ml2*m51*m52*xp2 A2*xpl-l/L*Blv*mll*ml2*m52*xp2 ; 
L34_31=-2/ (LA2) *Blv*mll*ml2*m52*xp2*xpl- 
1/ (LA3) *Blv*mll*ml2*m52*xp2*xplA2 ; 

L34_32=-l/ (2*L) *Blv*ml2 A2*m51*xp2-1/ (LA2) *Blv*ml2 A2*m51*xp2*xpl; 
L34_33=-l/ (2*LA3) *Blv*ml2A2*m51*xp2*xplA2- 
C1 v/ deltas at A2*kpsiA3*mll*ml2A2; 

L34_34=-Clv/ delta sat A2*kvA3*m21*m22 A2-Clv/ delta sat A2*kr A3*m31*m32 A2 ; 
L34_35=-Clv/ deltasat A2*kyA3*m41*m42 A2- 
Clv/ delta sat A2*kpsi A2*kv*m2 l*ml2 A2 ; 
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L34_36=-2*Clv/ delta sat A2*kpsiA2*kv*m22*mll*ml2- 
Civ/ delta sat A2*kpsiA2*kr*m31*ml2A2 ; 
L34_37=-2*Clv/deltasatA2*kpsiA2*kr*m32*mll*ml2- 
Civ/ delta sat A2*kpsiA2*ky*m4 l*ml2 A2 ; 
L34_38=-2*Clv/ delta sat A2*kpsi A2*ky*m42*ml l*ml2 ; 
L34_39=-2*Clv/ delta sat A2*kpsi*kvA2*m21*m22*ml2- 
Clv/ delta sat A2*kpsi*kvA2*m22 A2*ml 1 ; 
L34_40=-2*Clv/ delta sat A2*kpsi*kr A2*m31*m32*ml2- 


Clv/deltasatA2*kpsi* 
L34_41=-2*Clv/deltas 
Civ/ delta sat A2*kpsi* 
L34_42=-Clv/ delta sat 
2*Clv/ delta sat A2*kvA 
L34_43=-Clv/ delta sat 
2*Clv/ delta sat A2*kv 
L34_44=-Clv/ delta sat 
2*Clv/ delta sat A2*kv 
L34_45  =  -Clv/deltas 
2*Clv/ delta sat A2*kv* 
L34_46  =  -Clv/deltas 
2*Clv/ deltas at A2*kr 
L34_47  =  -Clv/deltas 
2*Clv/ deltas at A2*kr* 


L34_48  =  -2*Clv/delt 
L34_4  9  =  -2*Clv/delt 
L34_50  =  -2*Clv/delt 
L34_51  =  -2*Clv/ delt 
L34_52  =  -2*Clv/ delt 
L34_53  =  -2*Clv/ delt 
L34_54  =  -2*Clv/ delt 
L34_55  =  -2*Clv/ delt 
L34_56  =  -2*Clv/ delt 
L34_57  =  -2*Clv/ delt 
2*C1 v/ delt a sat A2*kv* 
L34_58  =  -2*Clv/ delt 


krA2*m32A2*mll ; 
at A2*kpsi*kyA2*m41*m42*ml2- 
kyA2*m42A2*mll ; 
A2*kvA2*kr*m31*m22A2- 
2*kr*m32*m21*m22 ; 

2*kvA2*ky*m41*m22A2- 
2*ky*m42*m21*m22 ; 

2*kv*krA2*m21*m32A2- 
kr/'2*m22*m31*m32 ; 
atA2*kv*kyA2*m21*m42A2- 
kyA2*m22*m41*m42 ; 
atA2*krA2*ky*m41*m32A2- 
2*ky*m42*m31*m32 ; 
atA2*kr*kyA2*m31*m42A2- 
kyA2*m32*m41*m42 ; 
a sat A2*kpsi*kv*kr*m31*m22*ml2 ; 
a sat A2*kpsi*kv*kr*m32*m2 1 *ml2 ; 
asatA2*kpsi*kv*kr*m32*m22*mll ; 
asatA2*kpsi*kv*ky*m41*m22*ml2 ; 
a sat A2*kpsi*kv*ky*m42*m2 1 *ml2 ; 
asatA2*kpsi*kv*ky*m42*m22*mll ; 
a sat A2*kpsi*kr*ky*m4 I*m32*ml2 ; 
asatA2*kpsi*kr*ky*m42*m31*ml2 ; 
asatA2*kpsi*kr*ky*m42*m32*mll ; 
a sat A2*kv*kr*ky*m41*m22*m32- 
kr*ky*m42*m21*m32; 
asatA2*kv*kr*ky*m42*m22*m31  ; 


L34  =  L34_l+L34_2+L34_3+L34_4+L34_5+L34_6+L34_7+L34_8+L34_9+L34_10 . . 
+L34_11+L34_12+L34_13+L34_14+L34_15+L34_16+L34_17+L34_18+L34_19+L34_20 
+L34_21+L34_22+L34_23+L34_24+L34_25+L34_26+L34_27+L34_28+L34_29+L34_30 
+L34_31+L34_32+L34_33+L34_34+L34_35+L34_36+L34_37+L34_38+L34_39+L34_40 
+L34_41+L34_42+L34_44+L34_45+L34_46+L34_47+L34_48+L34_49+L34_50+L34_51 
+L34_52+L34_53+L34_54+L34_55+L34_56+L34_57+L34_58; 


L4 1_1=1 / (6*LA3) *Blr*m41A3-l/ (6*LA3) *Blr*m81 A3-l/ (6*L) *Blr*m51"3*xp2; 
L4 1_2=-1 / (6*LA3) *Blr*m51A3*xp2A3-l/ (2*LA2) *Blr*mllA3*xplA2- 
1/ 6*Blr*mllA3; 

L41_3=-2/ (3*L) *Blr*mllA3*xpl- 

1/ (6*LA3) *Blr*mll"3*xplA3+l/ (2*LA3) *Blr*m41*m81A2; 
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L4 1 _ 4 =+ 1 / (2*LA3) *Blr*m41*m81A2-l/ (2*LA3) *Blr*xp2A2*m81*m51A2; 

L41_5  =  +1/ (2*LA3) *Blr*xp2 A2*m41*m51 A2-l/ (LA2) *Blr*m81*mllA2*xpl; 

L41_6  =  -1/ (2*LA3) *Blr*m81*mllA2*xplA2-l/ (2*L) *Blr*m81*mllA2; 

L4 1_7  =  -1/ (2*LA3) *Blr*xp2*m81A2*m51-l/ (2*LA2) *Blr*m81A2*mll; 

L4 1_8  =  -1/  (2*LA3)  *Blr*m81A2*mll*xpl-l/  (2*LA3)  *Blr*xp2*m41A2*m51; 

L41_9  =  -1/ (2*LA2) *Blr*m41A2*mll- 

1/ (2*LA3) *Blr*m41A2*mll*xpl+l/ (2*L) *Blr*m41*mllA2; 

L4 1_1 0  =  +1/ (LA2) *Blr*m41*mll A2*xpl+1/ (2*LA3) *Blr*m41*mllA2*xplA2; 

L4 1_1 1  =  +1/ (LA3) *Blr*xp2*m41*m81*m51+l/ (LA2) *Blr*m41*m81*mll; 

L41_12  =  +1/ (LA3) *Blr*m41*m81*mll*xpl-l/ (LA3) *Blr*m81*mll*m51*xp2*xpl; 
L41_13  =  -1/ (LA2) *Blr*m81*mll*m51*xp2+l/ (LA3) *Blr*m41*mll*m51*xp2*xpl; 

L4 1_1 4  =  +1/ (LA2) *Blr*m41*mll*m51*xp2-l/ (2*LA2) *Blr*mll*m51A2*xp2A2; 

L4 1_1 5  =  -1/ (2*LA3) *Blr*mll*m51 A2*xp2 A2*xpl-1/ (2*L) *Blr*mllA2*m51*xp2; 
L41_16  =  -1/ (LA2) *Blr*mllA2*m51*xp2*xpl- 
1/ (2*LA3) *Blr*mllA2*m51*xp2*xplA2 ; 

L4 1_1 7  =  -l/3*Clr/ deltas at A2*kpsiA3*mllA3-l/3*Clr/ deltasat A2*kvA3*m21 A3; 
L41_18  =  -l/3*Clr/ deltasat A 2 *krA 3 *m31A 3-1/ 3*Clr/ deltasat A2*kyA3*m41A3; 
L41_19  =  -Clr/deltasatA2*kpsiA2*kv*m21*mllA2~ 

Clr/ delta sat A2*kpsi A2*kr*m31*ml 1 A2 ; 

L41_20  =  -Clr/deltasatA2*kpsiA2*ky*m41*mllA2- 
Clr/ delta sat A2*kpsi*kvA2*m21 A2*ml 1 ; 

L41_21  =  -Clr/deltasat A2*kpsi*kr A2*m31 A2*mll- 
Clr/ delta sat A2*kpsi*kyA2*m4 1 A2*ml 1 ; 

L41_22  =  -Clr/ deltasat A2*kvA2*kr*m31*m21A2- 
Clr/ delta sat A2*kvA2*ky*m4 1 *m2 1 A2 ; 

L41_23  =  -Clr/ deltasat A2*kv*krA2*m21*m31A2- 
Clr/ delta sat A2*kv*kyA2*m21*m4 1 A2 ; 

L4 1_2  4  =  -Clr/ deltasat A2*krA2*ky*m41*m31A2- 
Clr/ delta sat A2*kr*kyA2*m31*m4 1 A2 ; 

L41_25  =  -2*Clr/deltasat A2*kpsi*kv*kr*m31 *m21 *ml 1 ; 

L41_26  =  -2*Clr/deltasat A2*kpsi*kv*ky*m4 l*m21*ml 1 ; 

L41_27  =  -2*Clr/deltasat A2*kpsi*kr*ky*m4 l*m31*ml 1 ; 

L41_28  =  -2*Clr/ deltasat A 2 *kv*kr*ky*m4 l*m2 l*m31 ; 


L4 1=L4 1_1+L4 1_2+L4 1_3+L4 1_4+L4 1_5+L4 1_6+L4 1_7  +L4 1_8+L4 1_9+L4 1_1 0 .  .  . 

+L4 1_1 1+L4 1_12+L4 1_1 3+L4 1_1 4+L4 1_1 5+L4 1_1 6+L4 1_1 7+L4 1_1 8+L4 1_1 9+L4 1_2  0 .  . 

+L41_21+L41_22+L41_23+L41_24+L41_25+L41_26+L41_27+L41_28; 

L42_l=l/ (6*LA3) *Blr*m42 A3-l/ (6*LA3) *Blr*m82 A3-l/ (6*L) *Blr*m52 A3*xp2  ; 
L42_2=-l/ (6*LA3) *Blr*m52 A3*xp2 A3-l/ (2*LA2) *Blr*ml2A3*xplA2-l/6*Blr*ml2A3 
L42_3=-2/ (3*L) *Blr*ml2 A3*xpl- 

1/ (6*LA3) *Blr*ml2A3*xplA3+l/ (2*LA3) *Blr*m42*m82 A2; 

L42_4=+l/ (2*LA3) *Blr*m42*m82 A2-l/ (2*LA3) *Blr*xp2 A2*m82*m52 A2  ; 

L42_5=+l/ (2*LA3) *Blr*xp2 A2*m42*m52 A2-l/ (LA2) *Blr*m82*ml2 A2*xpl; 

L42_6=-l/ (2*LA3) *Blr*m82*ml2 A2*xpl A2-l/ (2*L) *Blr*m82*ml2 A2  ; 

L42_7=-l/ (2*LA3) *Blr*xp2*m82 A2*m52-1/ (2*LA2) *Blr*m82 A2*ml2; 

L42_8=-l/ (2*LA3) *Blr*m82 A2*ml2*xpl-1/ (2*LA3) *Blr*xp2*m42 A2*m52; 

L42_9=-l/ (2*LA2 ) *Blr*m42 A2*ml2- 

1/ (2*LA3) *Blr*m42 A2*ml2*xpl+1/ (2*L) *Blr*m42*ml2 A2; 

L42_10=+l/ (LA2) *Blr*m42*ml2 A2*xpl+1/ (2*LA3) *Blr*m42*ml2 A2*xpl A2  ; 

L42_l 1  =  +1/ (LA3) *Blr*xp2*m42*m82*m52+l/ (LA2) *Blr*m42*m82*ml2; 
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L42_12  =  +1/ (LA3) *Blr*m42*m82*ml2*xpl-l/ (LA3) *Blr*m82*ml2*m52*xp2*xpl; 

L42_13  =  -1/ (LA2) *Blr*m82*ml2*m52*xp2+l/ (LA3) *Blr*m42*ml2*m52*xp2*xpl; 
L42_14=+l/ (LA2) *Blr*m42*ml2*m52*xp2-l/ (2*LA2) *Blr*ml2*m52 A2*xp2 A2  ; 
L42_15=-l/ (2*LA3) *Blr*ml2*m52 A2*xp2 A2*xpl-1/ (2*L) *Blr*ml2 A2*m52*xp2; 
L42_16=-l/ (LA2) *Blr*ml2 A2*m52*xp2*xpl-1/ (2*LA3) *Blr*ml2 A2*m52*xp2*xpl A2; 
L42_17=-l/ 3*Clr/ deltas at A2*kpsiA3*ml2A3-l/3*Clr/ deltasat A2*kvA3*m22 A3; 
L42_18  =-l /3*Clr /deltasat A2*krA3*m32A 3-1 /3*Clr /deltasat A2*kyA3*m4 2 A3  ; 

L42_19  =  -Clr/deltasat A2*kpsiA2*kv*m22*ml2 A2- 
Clr/ delta sat A2*kpsiA2*kr*m32*ml2 A2 ; 

L42_20  =  -Clr/deltasat A2*kpsiA2*ky*m42*ml2 A2- 
Clr/ delta sat A2*kpsi*kvA2*m22 A2*ml2 ; 

L42_21  =  -Clr/ deltasat A2*kpsi*kr A2*m32 A2*ml2- 
Clr/ delta sat A2*kpsi*kyA2*m42 A2*ml2 ; 

L42_22  =  -Clr/ deltasat A2*kvA2*kr*m32*m22A2- 
Clr/ delta sat A2*kvA2*ky*m42*m22 A2  ; 

L42_23  =  -Clr/ deltasat A2*kv*krA2*m22*m32A2- 
Clr/ delta sat A2*kv*kyA2*m22*m42 A2 ; 

L42_24  =-Clr/ deltasat A2*krA2*ky*m42*m32A2- 
Clr/ delta sat A2*kr*kyA2*m32*m42 A2 ; 

L42_25  =  -2*Clr/deltasat A2*kpsi*kv*kr*m32*m22*ml2; 

L42_26  =  -2*Clr/deltasatA2*kpsi*kv*ky*m42*m22*ml2; 

L42_27  =  -2*Clr/deltasat A2*kpsi*kr*ky*m42*m32*ml2; 

L42_28  =  -2*Clr/deltasat A2*kv*kr*ky*m42*m22*m32 ; 


L42=L42_1+L42_2+L42_3+L42_4+L42_5+L42_6+L42_7+L42_8+L42_9+L42_10 . . . 
+L42_11+L42_12+L42_13+L42_14+L42_15+L42_16+L42_17+L42_18+L42_19+L42_20 . . 
+L42_21+L42_22+L42_23+L42_24+L42_25+L42_26+L42_27+L42_28; 


L43_l=l/ (2*LA3) *Blr*m41 A2*m42-1/ (2*LA3) *Blr*m81 A2*m82- 
1/ (2*L) *Blr*m51 A2*m52*xp2 ; 

L43_2=-l/ (2*LA3) *Blr*m51 A2*m52*xp2 A3-3/ (2*LA2) *Blr*mllA2*ml2*xplA2  ; 

L43_3=-l/2*Blr*mllA2*ml2-2/L*Blr*mllA2*ml2*xpl- 

1/ (2*LA3) *Blr*mllA2*ml2*xplA3; 

L43_4=+l/ (LA3) *Blr*m41*m81*m82+l/ (2*LA3) *Blr*m42*m81 A2+l/ (LA3) *Blr*m41*m 
81*m82  ; 

L43_5=+l/ (2*LA3) *Blr*m42*m81A2-l/ (LA3) *Blr*xp2A2*m81*m51*m52; 

L43_6=-l/ (2*LA3) *Blr*xp2 A2*m82*m51 A2+l/ (LA3) *Blr*xp2A2*m41*m51*m52  ; 

L43_7=+l/ (2*LA3) *Blr*xp2 A2*m42*m51 A2-2/ (LA2) *Blr*m81*mll*ml2*xpl; 
L43_8=-l/ (LA3) *Blr*m81*mll*ml2*xplA2-l/L*Blr*m81*mll*ml2- 
1/ (LA2) *Blr*m82*mllA2*xpl  ; 

L43_9=-l/ (2*LA3) *Blr*m82*mllA2*xplA2-l/ (2*L) *Blr*m82*mll A2; 

L43_10=-l/ (2*LA3) *Blr*xp2*m81 A2*m52-1/ (LA3) *Blr*xp2*m81*m82*m51  ; 

L43_ll=-1/ (2*LA2 ) *Blr*m81 A2*ml2-1/ (2*LA3) *Blr*m81 A2*ml2*xpl- 
1/ (LA2) *Blr*m81*m82*mll; 

L43_12=-l/ (LA3) *Blr*m81*m82*mll*xpl-l/ (2*LA3) *Blr*xp2*m41A2*m52  ; 

L43_13=-l/ (LA3) *Blr*xp2*m41*m42*m51-l/ (2*LA2) *Blr*m41A2*ml2; 

L43_14=-l/ (2*LA3) *Blr*m41A2*ml2*xpl-l/ (LA2) *Blr*m41*m42*mll  ; 

L43_15=l/ (LA3) *Blr*m41*m42*mll*xpl+l/L*Blr*m41*mll*ml2+2/ (LA2) *Blr*m41*m 
ll*ml2*xpl; 

L43_16=+l/ (LA3) *Blr*m41*mll*ml2*xpl A2+l/ (2*L) *Blr*m42*mllA2  ; 
L43_17=+l/ (LA2) *Blr*m42*mllA2*xpl+l/ (2*LA3) *Blr*m42*mllA2*xplA2; 
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L43_18=+l/ (LA3) *Blr*xp2*m41*m81*m52+l/ (LA3) *Blr*xp2*m41*m82*m51  ; 
L43_19=+l/ (LA3) *Blr*xp2*m42*m81*m51+l/ (LA2) *Blr*m41*m81*ml2; 

L43_20=+l/ (LA3) *Blr*m41*m81*ml2*xpl+l/ (LA2) *Blr*m41*m82*mll  ; 

L43_21=+l/ (LA3) *Blr*m41*m82*mll*xpl+l/ (LA2) *Blr*m42*m81*mll; 

L43_22=+l/ (LA3) *Blr*m42*m81*mll*xpl-l/ (LA3) *Blr*m81*mll*m52*xp2*xpl  ; 

L43_23=-l/ (LA2) *Blr*m81*mll*m52*xp2-l/ (LA3) *Blr*m81*ml2*m51*xp2*xpl; 
L43_24=-l/ (LA2) *Blr*m81*ml2*m51*xp2-l/ (LA3) *Blr*m82*mll*m51*xp2*xpl  ; 
L43_25=-l/ (LA2) *Blr*m82*mll*m51*xp2+l/ (LA3) *Blr*m41*mll*m52*xp2*xpl; 
L43_26=+l/ (LA2) *Blr*m41*mll*m52*xp2+l/ (LA3) *Blr*m41*ml2*m51*xp2*xpl; 
L43_27=+l/ (LA2) *Blr*m41*ml2*m51*xp2+l/ (LA3) *Blr*m42*mll*m51*xp2*xpl; 
L43_28=+l/ (LA2) *Blr*m42*mll*m51*xp2-l/ (LA2) *Blr*mll*m51*m52*xp2 A2  ; 
L43_29=-l/ (LA3) *Blr*mll*m51*m52*xp2A2*xpl-l/ (2*LA2) *Blr*ml2*m51 A2*xp2 A2; 
L43_30=-l/ (2*LA3) *Blr*ml2*m51 A2*xp2 A2*xpl-1/ (2*L) *Blr*mll A2*m52*xp2  ; 

L43_31=-l/ (LA2) *Blr*mll A2*m52*xp2*xpl-1/ (2*LA3) *Blr*mllA2*m52*xp2*xplA2; 
L43_32=-l/L*Blr*mll*ml2*m51*xp2-2/ (LA2) *Blr*mll*ml2*m51*xp2*xpl  ; 

L43_33=-l/ (LA3) *Blr*mll*ml2*m51*xp2*xplA2- 
Clr/deltasatA2*kpsiA3*mllA2*ml2; 

L43_34=-Clr/ delta sat A2*kvA3*m21 A2*m22-Clr/ delta sat A2*kr A3*m31 A2*m32 ; 
L43_35=-Clr/ deltasat A2*kyA3*m41 A2*m42- 
2*Clr/ delta sat A2*kpsi A2*kv*m2 1 *ml l*ml2 ; 

L43_36=-Clr/ deltasat A2*kpsiA2*kv*m22*mll A2- 
2*Clr/ delta sat A2*kpsiA2*kr*m31*mll*ml2  ; 

L43_37=-Clr/ deltasat A2*kpsiA2*kr*m32*mll A2- 
2*Clr/ delta sat A2*kpsi A2*ky*m4 l*ml l*ml2 ; 

L43_38=-Clr/ delta sat A2*kpsiA2*ky*m42*mll A2- 
Clr/ delta sat A2*kpsi*kvA2*m21 A2*ml2  ; 

L43_39  =  -2*Clr/deltasat A2*kpsi*kvA2*m21*m22*mll- 
Clr/ delta sat A2*kpsi*kr A2*m31 A2*ml2 ; 

L43_40  =  -2*Clr/deltasatA2*kpsi*krA2*m31*m32*mll- 
Clr/ delta sat A2*kpsi*kyA2*m41 A2*ml2  ; 

L43_41  =-2*Clr/ delta sat A2*kpsi*kyA2*m41*m42*mll- 
2*Clr/ delta sat A2*kvA2*kr*m31*m21*m22  ; 

L43_42  =  -Clr/deltasat A2*kvA2*kr*m32*m21 A2- 
2*Clr/ delta sat A2*kvA2*ky*m41*m21*m22  ; 

L43_43  =  -Clr/ deltasat A2*kvA2*ky*m42*m21 A2- 
2*Clr/ delta sat A2*kv*kr A2*m21*m31*m32 ; 

L4  3_4  4  =  -Clr/ deltasat A2*kv*krA2*m22 *m31A2- 
2*Clr/ delta sat A2*kv*kyA2*m21*m4 l*m42 ; 

L43_45  =  -Clr/deltasatA2*kv*kyA2*m22*m41A2- 
2*Clr/ delta sat A2*kr A2*ky*m4 I*m31*m32 ; 

L43_46  =  -Clr/ deltasat A2*krA2*ky*m42 *m31A2- 
2*Clr/ delta sat A2*kr*kyA2*m31*m4 l*m42 ; 

L43_47  =  -Clr/ deltasat A2*kr*kyA2*m32*m41 A2- 
2*Clr/ delta sat A2*kpsi*kv*kr*m31*m21*ml2 ; 

L43_48  =  -2*Clr/deltasatA2*kpsi*kv*kr*m31*m22*mll  ; 

L43_49  =  -2*Clr/deltasat A2*kpsi*kv*kr*m32*m2 l*ml 1 ; 

L43_50  =  -2*Clr/deltasatA2*kpsi*kv*ky*m41*m21*ml2  ; 

L43_51  =  -2*Clr/deltasat A2*kpsi*kv*ky*m4 l*m22*ml 1 ; 

L43_52  =  -2*Clr/deltasat A2*kpsi*kv*ky*m42*m2 l*ml 1 ; 

L43_53  =  -2*Clr/deltasat A2*kpsi*kr*ky*m4 l*m31 *ml2 ; 

L43_54  =  -2*Clr/deltasatA2*kpsi*kr*ky*m41*m32*mll  ; 

L43_55  =  -2*Clr/deltasat A2*kpsi*kr*ky*m42*m31*ml 1 ; 

L43_56  =  -2*Clr/deltasatA2*kv*kr*ky*m41*m21*m32- 
2*Clr/ delta sat A2*kv*kr*ky*m4 I*m22*m31  ; 

L43_57  =  -2*Clr/ deltasat A 2 *kv*kr*ky*m4 2 *m2 l*m31 ; 
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L43=  L43_l+L43_2+L43_3+L43_4+L43_5+L43_6+L43_7+L43_8+L43_9+L43_10 . . . 

+L43_11+L43_12+L43_13+L43_14+L43_15+L43_16+L43_17+L43_18+L43_19+L43_20 . . 

+L43_21+L43_22+L43_23+L43_24+L43_25+L43_26+L43_27+L43_28+L43_29+L43_30 . . 

+L43_31+L43_32+L43_33+L43_34+L43_35+L43_36+L43_37+L43_38+L43_39+L43_40 . . 

+L43_41+L43_42+L43_44+L43_45+L43_46+L43_47+L43_48+L43_49+L43_50+L43_51 . . 
L43_52+L43_53+L43_54+L43_55+L43_56+L43_57; 


L4  4_1  =  1 / (2*LA3) *Blr*m41*m42A2-l/ (2*LA3) *Blr*m81*m82A2- 
1/ (2*L) *Blr*m51*m52 A2*xp2 ; 

L4  4_2=-l / (2*LA3) *Blr*m51*m52 A2*xp2 A3-3/ (2*LA2) *Blr*mll*ml2A2*xplA2; 

L44_3=-l/2*Blr*mll*ml2A2-2/L*Blr*mll*ml2A2*xpl- 

1/ (2*LA3) *Blr*mll*ml2 A2*xpl A3 ; 

L4  4_4=+l / (2*LA3) *Blr*m41*m82 A2+l/ (LA3) *Blr*m42*m81*m82+l/ (2*LA3) *Blr*m41 
*m82 A2 ; 

L44_5=+l/ (LA3)  *Blr*m42*m81*m82-l/ (2*LA3) *Blr*xp2A2*m81*m52A2; 

L44_6=-l/ (LA3) *Blr*xp2 A2*m82*m51*m52+1/ (2*LA3) *Blr*xp2A2*m41*m52A2; 
L44_7=+l/ (LA3) *Blr*xp2 A2*m42*m51*m52-1/ (LA2) *Blr*m81*ml2A2*xpl; 

L44_8=-l/ (2*LA3) *Blr*m81*ml2A2*xplA2-l/ (2*L) *Blr*m81*ml2A2; 

L44_9=-2/ (LA2) *Blr*m82*mll*ml2*xpl-l/ (LA3) *Blr*m82*mll*ml2*xplA2; 
L44_10=-l/L*Blr*m82*mll*ml2-l/ (LA3) *Blr*xp2*m81*m82*m52- 
1/ (2*LA3) *Blr*xp2*m82 A2*m51 ; 

L4  4_1 1=-1 / (LA2) *Blr*m81*m82*ml2-l/ (LA3) *Blr*m81*m82*ml2*xpl- 
1/ (2*LA2) *Blr*m82A2*mll  ; 

L44_12=— 1/ (2*LA3) *Blr*m82A2*mll*xpl-l/ (LA3) *Blr*xp2*m41*m42*m52; 
L44_13=-l/ (2*LA3) *Blr*xp2*m42 A2*m51-1/ (LA2) *Blr*m41*m42*ml2; 

L4  4_1 4=-l / (LA3) *Blr*m41*m42*ml2*xpl-l/ (2*LA2) *Blr*m42A2*mll; 

L44_15=-l/ (2*LA3) *Blr*m42 A2*mll*xpl+1/ (2*L) *Blr*m41*ml2 A2+ 

1/ (LA2) *Blr*m41*ml2A2*xpl ; 

L44_16=+l/ (2*LA3) *Blr*m41*ml2A2*xplA2+l/L*Blr*m42*mll*ml2; 

L4  4_1 7=+2 / (LA2) *Blr *m42 *ml 1 *ml2 *xpl+l / (LA3) *Blr*m42*mll*ml2*xplA2; 
L44_18=+l/ (LA3) *Blr*xp2*m41*m82*m52+l/ (LA3) *Blr*xp2*m42*m81*m52; 
L44_19=+l/ (LA3) *Blr*xp2*m42*m82*m51+l/ (LA2) *Blr*m41*m82*ml2; 

L44_20=+l/ (LA3) *Blr*m41*m82*ml2*xpl+l/ (LA2) *Blr*m42*m81*ml2; 

L44_21=+l/ (LA3) *Blr*m42*m81*ml2*xpl+l/ (LA2) *Blr*m42*m82*mll; 

L44_22=+l/ (LA3) *Blr*m42*m82*mll*xpl-l/ (LA3) *Blr*m81*ml2*m52*xp2*xpl; 
L44_23=-l/ (LA2) *Blr*m81*ml2*m52*xp2-l/ (LA3) *Blr*m82*mll*m52*xp2*xpl; 
L44_24=-l/ (LA2) *Blr *m82 *ml 1 *m52 *xp2-l / (LA3) *Blr*m82*ml2*m51*xp2*xpl; 

L4  4_2  5=-l / (LA2) *Blr*m82*ml2*m51*xp2+l/ (LA3) *Blr*m41*ml2*m52*xp2*xpl; 
L44_26=+l/ (LA2) *Blr*m41*ml2*m52*xp2+l/ (LA3) *Blr*m42*mll*m52*xp2*xpl; 
L44_27=+l/ (LA2) *Blr *m42 *ml 1 *m52 *xp2+l / (LA3) *Blr*m42*ml2*m51*xp2*xpl; 
L44_28=+l/ (LA2) *Blr*m42*ml2*m51*xp2-l/ (2*LA2) *Blr*mll*m52A2*xp2A2; 
L44_29=-l/ (2*LA3) *Blr*mll*m52A2*xp2A2*xpl-l/ (LA2) *Blr*ml2*m51*m52*xp2 A2; 
L44_30=-l/ (LA3) *Blr*ml2*m51*m52*xp2 A2*xpl-l/L*Blr*mll*ml2*m52*xp2 ; 
L44_31=-2/ (LA2) *Blr*mll*ml2*m52*xp2*xpl- 
1/ (LA3) *Blr*mll*ml2*m52*xp2*xpl A2 ; 

L44_32=-l/ (2*L) *Blr*ml2 A2*m51*xp2-1/ (LA2) *Blr*ml2 A2*m51*xp2*xpl; 
L44_33=-l/ (2*LA3) *Blr*ml2A2*m51*xp2*xplA2- 
Clr/ deltas at A2*kpsiA3*mll*ml2A2; 
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L44_34  =-Clr/ delta sat A2*kvA3*m21*m22 A2-Clr/ delta sat A2*krA3*m31*m32/v2 ; 

L44_35=-Clr/deltasatA2*kyA3*m41*m42A2- 

Clr/ delta sat A2*kpsi A2*kv*m2 l*ml2 A2 ; 

L44_36=-2*Clr/ delta sat A2*kpsiA2*kv*m22*mll*ml2- 
Clr/ delta sat A2*kpsiA2*kr*m31*ml2 A2 ; 
L44_37=-2*Clr/deltasatA2*kpsiA2*kr*m32*mll*ml2- 
Clr/ delta sat A2*kpsi/'2*ky*m4 l*ml2A2 ; 

L44_38=-2*Clr/ delta sat A2*kpsi A2*ky*m42*ml l*ml2 ; 

L44_39=-2*Clr/ delta sat A2*kpsi*kvA2*m21*m22*ml2- 
Clr/ delta sat A2*kpsi*kvA2*m22 A2*ml 1 ; 

L44_40  =  -2*Clr/deltasat A2*kpsi*kr A2*m31*m32*ml2- 
Clr/ delta sat A2*kpsi*kr A2*m32 A2*ml 1 ; 

L44_41=-2*Clr/ delta sat A2 *kpsi*kyA2 *m4 1 *m42 *ml2- 
Clr/ delta sat A2*kpsi*kyA2*m42 A2*ml 1 ; 

L44_42=-Clr/ deltasatA2*kvA2*kr*m31*m22A2- 
2*Clr/ delta sat A2*kvA2*kr*m32*m21*m22 ; 

L4  4_4  3=-Clr/  delta  sat  A2*kv/'2*ky*m41*m22/'2- 
2*Clr/ delta sat A2*kvA2*ky*m42*m2 l*m22 ; 

L44_44=-Clr/ deltasat A2 *kv*kr A2 *m2 1 *m32 A2- 
2*Clr/  deltasatA2*kv*kr''2*m22*m31*m32 ; 
L44_45=-Clr/deltasatA2*kv*kyA2*m21*m42A2- 
2*Clr/ delta sat A2*kv*kyA2*m22*m4 l*m42 ; 

L44_46=-Clr/ delta sat A2*kr A2*ky*m41*m32A2- 
2*Clr/ delta sat A2*krA2*ky*m42*m31*m32 ; 

L44_47=-Clr/ deltasat A2*kr*kyA2*m31*m42 A2- 
2*Clr/ delta sat A2*kr*kyA2*m32*m4 l*m42 ; 

L44_48=-2*Clr/ delta sat A2*kpsi*kv*kr*m31*m22*ml2 ; 

L44_49  =  -2*Clr/deltasat A2*kpsi*kv*kr*m32*m21*ml2; 

L44_50  =-2*Clr/ deltasatA2*kpsi*kv*kr*m32*m22*ml 1 ; 

L44_51  =  -2*Clr/deltasat A2*kpsi*kv*ky*m4 1 *m22*ml2 ; 

L44_52  =  -2*Clr/deltasatA2*kpsi*kv*ky*m42*m2 l*ml2 ; 

L44_53  =  -2*Clr/deltasat A2*kpsi*kv*ky*m42*m22*ml 1 ; 

L44_54  =-2*Clr/ deltasatA2*kpsi*kr*ky*m4 I*m32*ml2 ; 

L44_55  =  -2*Clr/deltasatA2*kpsi*kr*ky*m42*m31*ml2; 

L44_56  =  -2*Clr/deltasat A2*kpsi*kr*ky*m42*m32*ml 1 ; 

L44_57  =  -2*Clr/deltasatA2*kv*kr*ky*m41*m22*m32- 
2*Clr/ delta sat A2*kv*kr*ky*m42*m2 l*m32 ; 

L44_58  =-2*Clr/ deltasat/'2*kv*kr*ky*m42*m22*m31 ; 


L44  =  L4  4_1+L4  4_2+L4  4_3+L4  4_4+L4  4_5+L4  4_6+L4  4_7  +L4  4_8+L4  4_9  +  L4  4_1 0 .  . 

+L4  4_1 1+L4  4_12+L4  4_1 3+L4  4_1 4+L4  4_1 5+L4  4_1 6+L4  4_1 7  +L4  4_1 8+L4  4_1 9  +  L4  4_2  0 
+L4  4_2 1+L4  4_22+L4  4_2  3+L4  4_2  4+L4  4_2  5+L4  4_2  6+L4  4_2  7  +L4  4_2  8+L4  4_2  9  +  L4  4_3  0 
+L44_31+L44_32+L44_33+L44_34+L44_35+L44_36+L44_37+L44_38+L44_39+L44_40 
+L4  4_4 1+L4  4_42+L4  4_4  4+L4  4_4  5+L4  4_4  6+L4  4_4  7  +L4  4_4  8+L4  4_4  9+L4  4_50  +  L4  4_51 
+L44_52+L44_53+L44_54+L44_55+L44_56+L44_57+L44_58; 
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L51=l/6*mllA3-l/6*m51"3-l/2*m61*m51A21/2*m21*mllA2; 

L52=l/6*ml2A3-l/6*m52A3-l/2*m62*m52A2+l/2*m22*ml2A2; 

L53=l/2*mllA2*ml2-l/2*m51A2*m52-m61*m51*m52+m21*mll*ml2- 

l/2*m62*m51A2+l/2*m22*mllA2; 

L54=l/2*mll*ml2A2-l/2*m51*m52A2-l/2*m61*m52A2- 
m62*m51*m52+l/ 2*m21*ml2A2+m22*ml l*ml2 ; 


Rll 
R13 
R22 
R2  4 


N1 1 *L1 1  +  N12*L21 
N1 1 *L1 4  +  N12*L24 
N21*L13  +  N22*L23 
N21*L12  +  N22*L22 


+  N13*L31  +  N1 4  *L4 1  +  N15*L51; 
+  N13*L34  +  N1 4  *L4  4  +  N15*L54; 
+  N23*L33  +  N24*L43  +  N25*L53; 
+  N23*L32  +  N24*L42  +  N25*L52; 


Kcoef (index) = (1/8) * (3*R11+R13+R22+3*R24) ; 


end 


end 


end 


end 


%  Plot  (L,  crit  TC)  curve 
"6 

plot (L_v, Kcoef b ') , xlabel ( ' L ' ) , y label ('{ K} '), title ( ' Constant  T 
%end 


) , grid 
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